#Functions and packages

library(relaimpo)
## Loading required package: MASS
## Loading required package: boot
## Loading required package: survey
## Loading required package: grid
## Loading required package: Matrix
## Loading required package: survival
## 
## Attaching package: 'survival'
## The following object is masked from 'package:boot':
## 
##     aml
## 
## Attaching package: 'survey'
## The following object is masked from 'package:graphics':
## 
##     dotchart
## Loading required package: mitools
## This is the global version of package relaimpo.
## If you are a non-US user, a version with the interesting additional metric pmvd is available
## from Ulrike Groempings web site at prof.beuth-hochschule.de/groemping.
library(ggplot2)
library(ggtext)
library(plyr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(reshape2)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library(cowplot)
library(lme4)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:boot':
## 
##     logit
library(doBy)
## 
## Attaching package: 'doBy'
## The following object is masked from 'package:dplyr':
## 
##     order_by
library(bipartite)
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## 
## Attaching package: 'lattice'
## The following object is masked from 'package:boot':
## 
##     melanoma
## This is vegan 2.5-7
## 
## Attaching package: 'vegan'
## The following object is masked from 'package:survey':
## 
##     calibrate
## Loading required package: sna
## Loading required package: statnet.common
## 
## Attaching package: 'statnet.common'
## The following objects are masked from 'package:base':
## 
##     attr, order
## Loading required package: network
## 
## 'network' 1.17.1 (2021-06-12), part of the Statnet Project
## * 'news(package="network")' for changes since last version
## * 'citation("network")' for citation information
## * 'https://statnet.org' for help, support, and other information
## 
## Attaching package: 'network'
## The following object is masked from 'package:plyr':
## 
##     is.discrete
## sna: Tools for Social Network Analysis
## Version 2.6 created on 2020-10-5.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
##  For citation information, type citation("sna").
##  Type help(package="sna") to get started.
##  This is bipartite 2.16.
##  For latest changes see versionlog in ?"bipartite-package". For citation see: citation("bipartite").
##  Have a nice time plotting and analysing two-mode networks.
## 
## Attaching package: 'bipartite'
## The following object is masked from 'package:vegan':
## 
##     nullmodel
## The following object is masked from 'package:plyr':
## 
##     empty
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ tibble  3.1.4     ✓ purrr   0.3.4
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   2.0.1     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::arrange()     masks plyr::arrange()
## x gridExtra::combine() masks dplyr::combine()
## x purrr::compact()     masks plyr::compact()
## x dplyr::count()       masks plyr::count()
## x tidyr::expand()      masks Matrix::expand()
## x dplyr::failwith()    masks plyr::failwith()
## x dplyr::filter()      masks stats::filter()
## x dplyr::id()          masks plyr::id()
## x dplyr::lag()         masks stats::lag()
## x dplyr::mutate()      masks plyr::mutate()
## x doBy::order_by()     masks dplyr::order_by()
## x tidyr::pack()        masks Matrix::pack()
## x car::recode()        masks dplyr::recode()
## x dplyr::rename()      masks plyr::rename()
## x dplyr::select()      masks MASS::select()
## x purrr::some()        masks car::some()
## x dplyr::summarise()   masks plyr::summarise()
## x dplyr::summarize()   masks plyr::summarize()
## x tidyr::unpack()      masks Matrix::unpack()
library(ggforce)
library(ggpubr)
## 
## Attaching package: 'ggpubr'
## The following object is masked from 'package:cowplot':
## 
##     get_legend
## The following object is masked from 'package:plyr':
## 
##     mutate
library(ggeffects)
## 
## Attaching package: 'ggeffects'
## The following object is masked from 'package:cowplot':
## 
##     get_title
library(viridis)
## Loading required package: viridisLite
library(viridisLite)
library(phia)
library(jtools)
library(EcoSimR)
library(cooccur)
library(fitdistrplus)

  
theme_ines<-theme(axis.text = element_text(size=14), axis.title = element_text(size=14, face="bold"), legend.text = element_text(size=12), strip.text = element_text(size=14), plot.title = element_text(size=14, face="bold"), panel.grid=element_line(colour="white"), panel.background = element_rect(fill="white") , axis.line = element_line(size = 0.5, linetype = "solid",
                                   colour = "black"), strip.background = element_rect(fill="white"))

save_plot<-function(dir, width=15, height=10, ...){
  ggsave(dir, width = width, height = height, units = c("cm"))
}

1 - Importing data

Importing coexistence data and the leaf age

df<-read.csv("./Data/Coexistence_data_final_corr.csv", header=TRUE, row.names = NULL, as.is=c(2))
dir.create("./Plots")
## Warning in dir.create("./Plots"): './Plots' already exists
leaf_conversion<-read.table("./Data/leaf_conversion.txt", header=TRUE, row.names = NULL, as.is=c(1,4))

df$leaf_age<-sapply(c(1:length(df[,1])), function(x){
  leaf_conversion$Age[which(as.character(leaf_conversion$Treatment)==as.character(df$Treatment[x]) & leaf_conversion$Box==df$Box[x] & leaf_conversion$Leaf==df$Leaf[x])]
})

df$Density<-mapvalues(df$Treatment, as.character(levels(as.factor(df$Treatment))), c("10:10","10:10", "10:10", "19:1", "19:1", "1:19","1:19", "20:0", "0:20"))

df$Timing<-mapvalues(df$Treatment, as.character(levels(as.factor(df$Treatment))), c("same time", "Tu first", "Te first", "same time", "Te first", "same time", "Tu first", "single", "single"))

str(df)
## 'data.frame':    3214 obs. of  11 variables:
##  $ Box      : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Treatment: chr  "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
##  $ Block    : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Leaf     : int  2 2 2 2 2 2 3 3 3 3 ...
##  $ Leaflet  : int  1 2 3 1 2 3 1 2 3 4 ...
##  $ Direction: Factor w/ 2 levels "Down","Up": 1 1 1 2 2 2 1 1 1 1 ...
##  $ Te       : int  4 2 1 6 6 3 0 0 4 0 ...
##  $ Tu       : int  1 2 5 1 12 6 1 2 0 2 ...
##  $ leaf_age : chr  "old" "old" "old" "old" ...
##  $ Density  : chr  "10:10" "10:10" "10:10" "10:10" ...
##  $ Timing   : chr  "same time" "same time" "same time" "same time" ...

Calculating proportions and ratios

df$ratio<-sapply(c(1:length(df[,1])), function(x){
  #print(x)
  if(df$Treatment[x]=="20Te" | df$Treatment[x]=="20Tu"){
  a<-NA
  }else{
  if(df$Te[x]==0){
    te<-1
  }else{
    te<-df$Te[x]
  }
  
  if(df$Tu[x]==0){
    tu<-1
  }else{
    tu<-df$Tu[x]
  }
    
     a<-te/tu
  }

  a
})


df$Total<-sapply(c(1:length(df[,1])), function(x){
  if(df$Treatment[x]=="20Te"){
  a<-df$Te[x]
  }else if(df$Treatment[x]=="20Tu"){
    a<-df$Tu[x]
  }else{
     a<-df$Te[x]+df$Tu[x]
  }

  a
})

df$proportionTe<-sapply(c(1:length(df[,1])), function(x){
  if(df$Treatment[x]=="20Tu"){
  a<-0
  }else if(df$Te[x]==0 |df$Total[x]==0 ){
    a<-0
    } else{
    te<-df$Te[x]
    a<-te/df$Total[x]
  }

  a
})

df$proportionTu<-sapply(c(1:length(df[,1])), function(x){
  if(df$Treatment[x]=="20Te"){
  a<-0
  }else if(df$Tu[x]==0 |df$Total[x]==0 ){
    a<-0
    } else{
    tu<-df$Tu[x]
    a<-tu/df$Total[x]
  }

  a
})

Setting up data frame with data per box

For the statistical analysis we will analyse data per box first. Here we are summing the number of individuals of each species per box - df_box To test a different distribution per leaf we will also create a data frame with number of individuals per leaf and per direction - df_box_leaf_dir

df_box_leaf_dir<-df %>%
  group_by(Treatment, Density, Timing , Block, Box, Direction,Leaf, leaf_age) %>%
     summarise(Te=sum(Te, na.rm=TRUE), Tu=sum(Tu, na.rm=TRUE)) %>%
  as.data.frame()
## `summarise()` has grouped output by 'Treatment', 'Density', 'Timing', 'Block', 'Box', 'Direction', 'Leaf'. You can override using the `.groups` argument.
df_box<-df %>%
  group_by(Treatment,Density, Timing ,  Block, Box) %>%
  summarise(Te=sum(Te, na.rm=TRUE), Tu=sum(Tu, na.rm=TRUE)) %>%
  as.data.frame()
## `summarise()` has grouped output by 'Treatment', 'Density', 'Timing', 'Block'. You can override using the `.groups` argument.
head(df_box_leaf_dir)
## Lets add our usual proportion information

df_box$ratio<-sapply(c(1:length(df_box[,1])), function(x){
  #print(x)
  if(df_box$Treatment[x]=="20Te" | df_box$Treatment[x]=="20Tu"){
  a<-NA
  }else{
  if(df_box$Te[x]==0){
    te<-1
  }else{
    te<-df_box$Te[x]
  }
  
  if(df_box$Tu[x]==0){
    tu<-1
  }else{
    tu<-df_box$Tu[x]
  }
    
     a<-te/tu
  }

  a
})

df_box_leaf_dir$ratio<-sapply(c(1:length(df_box_leaf_dir[,1])), function(x){
  #print(x)
  if(df_box_leaf_dir$Treatment[x]=="20Te" | df_box_leaf_dir$Treatment[x]=="20Tu"){
  a<-NA
  }else{
  if(df_box_leaf_dir$Te[x]==0){
    te<-1
  }else{
    te<-df_box_leaf_dir$Te[x]
  }
  
  if(df_box_leaf_dir$Tu[x]==0){
    tu<-1
  }else{
    tu<-df_box_leaf_dir$Tu[x]
  }
    
     a<-te/tu
  }

  a
})



df_box_leaf_dir$Total<-sapply(c(1:length(df_box_leaf_dir[,1])), function(x){
  if(df_box_leaf_dir$Treatment[x]=="20Te"){
  a<-df_box_leaf_dir$Te[x]
  }else if(df_box_leaf_dir$Treatment[x]=="20Tu"){
    a<-df_box_leaf_dir$Tu[x]
  }else{
     a<-df_box_leaf_dir$Te[x]+df_box_leaf_dir$Tu[x]
  }

  a
})


df_box_leaf_dir$proportionTe<-sapply(c(1:length(df_box_leaf_dir[,1])), function(x){
  if(df_box_leaf_dir$Treatment[x]=="20Tu"){
  a<-0
  }else if(df_box_leaf_dir$Te[x]==0 |df_box_leaf_dir$Total[x]==0 ){
    a<-0
    } else{
    te<-df_box_leaf_dir$Te[x]
    a<-te/df_box_leaf_dir$Total[x]
  }

  a
})

df_box_leaf_dir$proportionTu<-sapply(c(1:length(df_box_leaf_dir[,1])), function(x){
  if(df_box_leaf_dir$Treatment[x]=="20Te"){
  a<-0
  }else if(df_box_leaf_dir$Tu[x]==0 |df_box_leaf_dir$Total[x]==0 ){
    a<-0
    } else{
    tu<-df_box_leaf_dir$Tu[x]
    a<-tu/df_box_leaf_dir$Total[x]
  }

  a
})


df_box$Total<-sapply(c(1:length(df_box[,1])), function(x){
  if(df_box$Treatment[x]=="20Te"){
  a<-df_box$Te[x]
  }else if(df_box$Treatment[x]=="20Tu"){
    a<-df_box$Tu[x]
  }else{
     a<-df_box$Te[x]+df_box$Tu[x]
  }

  a
})


df_box$proportionTe<-sapply(c(1:length(df_box[,1])), function(x){
  if(df_box$Treatment[x]=="20Tu"){
  a<-0
  }else if(df_box$Te[x]==0 |df_box$Total[x]==0 ){
    a<-0
    } else{
    te<-df_box$Te[x]
    a<-te/df_box$Total[x]
  }

  a
})

df_box$proportionTu<-sapply(c(1:length(df_box[,1])), function(x){
  if(df_box$Treatment[x]=="20Te"){
  a<-0
  }else if(df_box$Tu[x]==0 |df_box$Total[x]==0 ){
    a<-0
    } else{
    tu<-df_box$Tu[x]
    a<-tu/df_box$Total[x]
  }

  a
})

Checking distribution of data

The variables that we are interested to know the distribution are: Te and Tu (because of the binomial analyses)

hist(df$Te)

hist(df$Tu)

descdist(df$Te[-which(is.na(df$Te/df$Total))], discrete = FALSE, boot=500)

## summary statistics
## ------
## min:  0   max:  149 
## median:  4 
## mean:  10.08053 
## estimated sd:  15.99749 
## estimated skewness:  3.468693 
## estimated kurtosis:  19.60968
descdist(df$Tu[-which(is.na(df$Tu/df$Total))], discrete = FALSE, boot=500)

## summary statistics
## ------
## min:  0   max:  125 
## median:  2 
## mean:  4.420565 
## estimated sd:  8.373305 
## estimated skewness:  5.241611 
## estimated kurtosis:  47.52895

Clearly not a normal distribution so we will use the binomial family in glmer to analyse data

Statistical analyses

2 - Does the timing and initial density change the probability of coexistence?

2.1 - Whole data set

Figure 1

df_box$Density2<-factor(df_box$Density, levels=c("1:19","10:10","19:1"))



  ggplot(subset(df_box, Treatment!="20Te" & Treatment!="20Tu"), aes(x=Density2, y=proportionTe, fill=interaction(Density2,Timing)))+
    facet_wrap("Timing", ncol=3, scales = "free_x")+
    geom_hline(yintercept = 0.5)+
    geom_boxplot()+
    theme_ines+
    scale_fill_manual(values = c("#e0f3f8", "#91bfdb", "#4575b4","#d73027", "#fc8d59", "#fee090", "#ffffbf"), name="Treatment", guide="none", drop=TRUE)+
    ylab(c("Proportion of T. evansi"))+
    xlab(c("Initial T. evansi: T. urticae density"))+
    #scale_x_discrete(labels=c())+
    theme(legend.title = element_text(size=18, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75))+
    stat_summary(fun.y=median, geom="smooth", aes(group=Timing, colour=Timing), lwd=1)+
    scale_colour_manual(values = c( "#4575b4","#d73027", "#fee090"), name="Treatment", guide="none")
## Warning: `fun.y` is deprecated. Use `fun` instead.

save_plot("./Plots/Figure1.pdf", width=15, height=10)

Stats

Random part of the model

First we need to remove the single species regimes, to test for the effects Since here we only have as random variable Block, our model will include block as random factor Our unit of analysis will be the box

df$Block2<-as.factor(df$Block)
df$Box2<-as.factor(df$Box)
str(df)
## 'data.frame':    3214 obs. of  17 variables:
##  $ Box         : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Treatment   : chr  "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
##  $ Block       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Leaf        : int  2 2 2 2 2 2 3 3 3 3 ...
##  $ Leaflet     : int  1 2 3 1 2 3 1 2 3 4 ...
##  $ Direction   : Factor w/ 2 levels "Down","Up": 1 1 1 2 2 2 1 1 1 1 ...
##  $ Te          : int  4 2 1 6 6 3 0 0 4 0 ...
##  $ Tu          : int  1 2 5 1 12 6 1 2 0 2 ...
##  $ leaf_age    : chr  "old" "old" "old" "old" ...
##  $ Density     : chr  "10:10" "10:10" "10:10" "10:10" ...
##  $ Timing      : chr  "same time" "same time" "same time" "same time" ...
##  $ ratio       : num  4 1 0.2 6 0.5 0.5 1 0.5 4 0.5 ...
##  $ Total       : int  5 4 6 7 18 9 1 2 4 2 ...
##  $ proportionTe: num  0.8 0.5 0.167 0.857 0.333 ...
##  $ proportionTu: num  0.2 0.5 0.833 0.143 0.667 ...
##  $ Block2      : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Box2        : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
df_wsingle<-subset(df, Treatment!="20Te" & Treatment!="20Tu")
df_wsingle<-droplevels(df_wsingle)
str(df_wsingle)
## 'data.frame':    2500 obs. of  17 variables:
##  $ Box         : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Treatment   : chr  "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
##  $ Block       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Leaf        : int  2 2 2 2 2 2 3 3 3 3 ...
##  $ Leaflet     : int  1 2 3 1 2 3 1 2 3 4 ...
##  $ Direction   : Factor w/ 2 levels "Down","Up": 1 1 1 2 2 2 1 1 1 1 ...
##  $ Te          : int  4 2 1 6 6 3 0 0 4 0 ...
##  $ Tu          : int  1 2 5 1 12 6 1 2 0 2 ...
##  $ leaf_age    : chr  "old" "old" "old" "old" ...
##  $ Density     : chr  "10:10" "10:10" "10:10" "10:10" ...
##  $ Timing      : chr  "same time" "same time" "same time" "same time" ...
##  $ ratio       : num  4 1 0.2 6 0.5 0.5 1 0.5 4 0.5 ...
##  $ Total       : int  5 4 6 7 18 9 1 2 4 2 ...
##  $ proportionTe: num  0.8 0.5 0.167 0.857 0.333 ...
##  $ proportionTu: num  0.2 0.5 0.833 0.143 0.667 ...
##  $ Block2      : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Box2        : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
levels(as.factor(df_wsingle$Treatment))
## [1] "10Te10Tu"    "10Te10Tu48h" "10Te48h10Tu" "19Te1Tu"     "19Te48h1Tu" 
## [6] "1Te19Tu"     "1Te19Tu48h"
r1<-glmer(cbind(Te,Total)~1+ (1|Block2), data=df_wsingle, family = binomial(link="logit"))

summary(r1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(Te, Total) ~ 1 + (1 | Block2)
##    Data: df_wsingle
## 
##      AIC      BIC   logLik deviance df.resid 
##  11424.3  11436.0  -5710.2  11420.3     2498 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -9.0545 -0.8376  0.0000  0.5633  3.0748 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev.
##  Block2 (Intercept) 0.0002307 0.01519 
## Number of obs: 2500, groups:  Block2, 2
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.3659     0.0143  -25.58   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Complete model with Block as random factor
b1<-glmer(cbind(Te,Total)~ Density*Timing+(1|Block2), data=df_wsingle, family = binomial(link="logit"))
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
## boundary (singular) fit: see ?isSingular
b2<-glmer(cbind(Te,Total)~ Density+Timing+(1|Block2), data=df_wsingle, family = binomial(link="logit"))
## boundary (singular) fit: see ?isSingular
b3<-glmer(cbind(Te,Total)~ Timing+(1|Block2), data=df_wsingle, family = binomial(link="logit"))
b4<-glmer(cbind(Te,Total)~ Density+(1|Block2), data=df_wsingle, family = binomial(link="logit"))
## boundary (singular) fit: see ?isSingular
b5<-glmer(cbind(Te,Total)~ Treatment+(1|Block2), data=df_wsingle, family = binomial(link="logit"))
## boundary (singular) fit: see ?isSingular
anova(b1, b2, b3, b4, b5)

The lowest AIC is from models b1 and b5, which correspond to the complete model, so lets take a look at the summary of the model

summary(b1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(Te, Total) ~ Density * Timing + (1 | Block2)
##    Data: df_wsingle
## 
##      AIC      BIC   logLik deviance df.resid 
##   7602.9   7649.4  -3793.4   7586.9     2492 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.3247 -0.3826  0.0000  0.1353  9.5786 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev. 
##  Block2 (Intercept) 2.142e-17 4.628e-09
## Number of obs: 2500, groups:  Block2, 2
## 
## Fixed effects:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 -1.30669    0.05790 -22.567  < 2e-16 ***
## Density10:10                 0.81807    0.06352  12.880  < 2e-16 ***
## Density19:1                  1.28768    0.06087  21.155  < 2e-16 ***
## TimingTe first              -0.01084    0.02572  -0.421    0.673    
## TimingTu first              -1.71343    0.09913 -17.285  < 2e-16 ***
## Density10:10:TimingTe first  0.30846    0.04311   7.155 8.37e-13 ***
## Density10:10:TimingTu first  1.35317    0.10696  12.651  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                            (Intr) D10:10 Dn19:1 TimingTefirst TimingTufirst
## Densty10:10                -0.912                                          
## Density19:1                -0.951  0.867                                   
## TimingTefirst               0.000  0.000 -0.225                            
## TimingTufirst              -0.584  0.533  0.556  0.000                     
## Density10:10:TimingTefirst  0.000 -0.249  0.134 -0.597         0.000       
## Density10:10:TimingTufirst  0.541 -0.594 -0.515  0.000        -0.927       
##                            Density10:10:TimingTefirst
## Densty10:10                                          
## Density19:1                                          
## TimingTefirst                                        
## TimingTufirst                                        
## Density10:10:TimingTefirst                           
## Density10:10:TimingTufirst  0.148                    
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficients
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
summary(b5)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(Te, Total) ~ Treatment + (1 | Block2)
##    Data: df_wsingle
## 
##      AIC      BIC   logLik deviance df.resid 
##   7602.9   7649.4  -3793.4   7586.9     2492 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.3247 -0.3826  0.0000  0.1353  9.5786 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Block2 (Intercept) 0        0       
## Number of obs: 2500, groups:  Block2, 2
## 
## Fixed effects:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -0.48862    0.02610 -18.719   <2e-16 ***
## Treatment10Te10Tu48h -0.36026    0.04017  -8.967   <2e-16 ***
## Treatment10Te48h10Tu  0.29763    0.03460   8.602   <2e-16 ***
## Treatment19Te1Tu      0.46961    0.03215  14.608   <2e-16 ***
## Treatment19Te48h1Tu   0.45877    0.03148  14.575   <2e-16 ***
## Treatment1Te19Tu     -0.81807    0.06352 -12.880   <2e-16 ***
## Treatment1Te19Tu48h  -2.53150    0.08459 -29.927   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) T10T10 T10T48 T19T1T T19T48 Tr1T19T
## Trt10T10T48 -0.650                                    
## Trt10T4810T -0.754  0.490                             
## Trtmnt19T1T -0.812  0.528  0.613                      
## Trtm19T481T -0.829  0.539  0.626  0.673               
## Trtmnt1T19T -0.411  0.267  0.310  0.334  0.341        
## Trtm1T19T48 -0.309  0.200  0.233  0.251  0.256  0.127 
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
Anova(b1, type=3)
Anova(b5, type=3)

There is a clear effect of treatment in the proportion of Te females. To test the impact of orde of arrival and initial frequency we can perform contrasts using the phia package.

Contrasts

Using a single factor

b5_2<-glmer(cbind(Te,Total)~ Treatment+(1|Block2), data=df_wsingle, family = binomial(link="logit"))
## boundary (singular) fit: see ?isSingular
summary(b5_2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(Te, Total) ~ Treatment + (1 | Block2)
##    Data: df_wsingle
## 
##      AIC      BIC   logLik deviance df.resid 
##   7602.9   7649.4  -3793.4   7586.9     2492 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.3247 -0.3826  0.0000  0.1353  9.5786 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Block2 (Intercept) 0        0       
## Number of obs: 2500, groups:  Block2, 2
## 
## Fixed effects:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -0.48862    0.02610 -18.719   <2e-16 ***
## Treatment10Te10Tu48h -0.36026    0.04017  -8.967   <2e-16 ***
## Treatment10Te48h10Tu  0.29763    0.03460   8.602   <2e-16 ***
## Treatment19Te1Tu      0.46961    0.03215  14.608   <2e-16 ***
## Treatment19Te48h1Tu   0.45877    0.03148  14.575   <2e-16 ***
## Treatment1Te19Tu     -0.81807    0.06352 -12.880   <2e-16 ***
## Treatment1Te19Tu48h  -2.53150    0.08459 -29.927   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) T10T10 T10T48 T19T1T T19T48 Tr1T19T
## Trt10T10T48 -0.650                                    
## Trt10T4810T -0.754  0.490                             
## Trtmnt19T1T -0.812  0.528  0.613                      
## Trtm19T481T -0.829  0.539  0.626  0.673               
## Trtmnt1T19T -0.411  0.267  0.310  0.334  0.341        
## Trtm1T19T48 -0.309  0.200  0.233  0.251  0.256  0.127 
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
levels(as.factor(df_wsingle$Treatment))
## [1] "10Te10Tu"    "10Te10Tu48h" "10Te48h10Tu" "19Te1Tu"     "19Te48h1Tu" 
## [6] "1Te19Tu"     "1Te19Tu48h"
test_int2<-testInteractions(b5_2, pairwise="Treatment", adjustment = "fdr")
test_int2
#Now we just select the ones we want to compare
#effect of timing
test_int2[c(1, 2,7,16, 21 ),]
#effect of density
test_int2[c(3,5,11, 13,17 ),]
# Testing overall order of arrival and overall density

mat_overall<-matrix(ncol=4, nrow=7)
colnames(mat_overall)<-c("Order of arrival Te first", "Order of arrival Tu first", "Te higher initially", "Tu higher initially")
rownames(mat_overall)<-c("10Te10Tu","10Te10Tu48h", "10Te48h10Tu", "19Te1Tu", "19Te48h1Tu", "1Te19Tu", "1Te19Tu48h" )

mat_overall[1,1]<- 1
mat_overall[2,1]<- 0
mat_overall[3,1]<- -1
mat_overall[4,1]<- 1
mat_overall[5,1]<- -1
mat_overall[6,1]<- 0
mat_overall[7,1]<- 0

mat_overall[1,2]<- 1
mat_overall[2,2]<- -1
mat_overall[3,2]<- 0
mat_overall[4,2]<- 0
mat_overall[5,2]<- 0
mat_overall[6,2]<- 1
mat_overall[7,2]<- -1


mat_overall[1,3]<- 1
mat_overall[2,3]<- 1
mat_overall[3,3]<- 1
mat_overall[4,3]<- -1
mat_overall[5,3]<- -1
mat_overall[6,3]<- 0
mat_overall[7,3]<- 0

mat_overall[1,4]<- 1
mat_overall[2,4]<- 1
mat_overall[3,4]<- 1
mat_overall[4,4]<- 0
mat_overall[5,4]<- 0
mat_overall[6,4]<- -1
mat_overall[7,4]<- -1

mat_overall
##             Order of arrival Te first Order of arrival Tu first
## 10Te10Tu                            1                         1
## 10Te10Tu48h                         0                        -1
## 10Te48h10Tu                        -1                         0
## 19Te1Tu                             1                         0
## 19Te48h1Tu                         -1                         0
## 1Te19Tu                             0                         1
## 1Te19Tu48h                          0                        -1
##             Te higher initially Tu higher initially
## 10Te10Tu                      1                   1
## 10Te10Tu48h                   1                   1
## 10Te48h10Tu                   1                   1
## 19Te1Tu                      -1                   0
## 19Te48h1Tu                   -1                   0
## 1Te19Tu                       0                  -1
## 1Te19Tu48h                    0                  -1
test_int3<-testInteractions(b5_2, custom=list(Treatment=mat_overall), adjustment = "fdr")
test_int3

2.2 - Test the importance of order of leaves

To test this we will run the model for the subset of the data with leaves 2 -4 first or 3 - 5 first

str(df_wsingle)
## 'data.frame':    2500 obs. of  17 variables:
##  $ Box         : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Treatment   : chr  "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
##  $ Block       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Leaf        : int  2 2 2 2 2 2 3 3 3 3 ...
##  $ Leaflet     : int  1 2 3 1 2 3 1 2 3 4 ...
##  $ Direction   : Factor w/ 2 levels "Down","Up": 1 1 1 2 2 2 1 1 1 1 ...
##  $ Te          : int  4 2 1 6 6 3 0 0 4 0 ...
##  $ Tu          : int  1 2 5 1 12 6 1 2 0 2 ...
##  $ leaf_age    : chr  "old" "old" "old" "old" ...
##  $ Density     : chr  "10:10" "10:10" "10:10" "10:10" ...
##  $ Timing      : chr  "same time" "same time" "same time" "same time" ...
##  $ ratio       : num  4 1 0.2 6 0.5 0.5 1 0.5 4 0.5 ...
##  $ Total       : int  5 4 6 7 18 9 1 2 4 2 ...
##  $ proportionTe: num  0.8 0.5 0.167 0.857 0.333 ...
##  $ proportionTu: num  0.2 0.5 0.833 0.143 0.667 ...
##  $ Block2      : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Box2        : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
pair2_4_old<-subset(df_wsingle, (Leaf==2 & leaf_age=="old") | (Leaf==4 & leaf_age=="old") | (Leaf==3 & leaf_age=="new") | (Leaf==5 & leaf_age=="new"))

pair2_4_new<-subset(df_wsingle, (Leaf==3 & leaf_age=="old") | (Leaf==5 & leaf_age=="old") | (Leaf==2 & leaf_age=="new") | (Leaf==4 & leaf_age=="new"))

g24Old<-glmer(cbind(Te,Total)~ Treatment+(1|Block2), data=pair2_4_old, family = binomial(link="logit"))
## boundary (singular) fit: see ?isSingular
summary(g24Old)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(Te, Total) ~ Treatment + (1 | Block2)
##    Data: pair2_4_old
## 
##      AIC      BIC   logLik deviance df.resid 
##   3768.4   3809.5  -1876.2   3752.4     1242 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4200 -0.4321  0.0000  0.1346  9.7348 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Block2 (Intercept) 0        0       
## Number of obs: 1250, groups:  Block2, 2
## 
## Fixed effects:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -0.52429    0.03926 -13.355  < 2e-16 ***
## Treatment10Te10Tu48h -0.11894    0.05602  -2.123   0.0337 *  
## Treatment10Te48h10Tu  0.32034    0.05050   6.343 2.25e-10 ***
## Treatment19Te1Tu      0.50774    0.04666  10.881  < 2e-16 ***
## Treatment19Te48h1Tu   0.49115    0.04790  10.254  < 2e-16 ***
## Treatment1Te19Tu     -0.86836    0.09743  -8.913  < 2e-16 ***
## Treatment1Te19Tu48h  -2.52471    0.11774 -21.442  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) T10T10 T10T48 T19T1T T19T48 Tr1T19T
## Trt10T10T48 -0.701                                    
## Trt10T4810T -0.777  0.545                             
## Trtmnt19T1T -0.841  0.590  0.654                      
## Trtm19T481T -0.820  0.574  0.637  0.689               
## Trtmnt1T19T -0.403  0.282  0.313  0.339  0.330        
## Trtm1T19T48 -0.333  0.234  0.259  0.280  0.273  0.134 
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
Anova(g24Old, type=3)
g24New<-glmer(cbind(Te,Total)~ Treatment+(1|Block2), data=pair2_4_new, family = binomial(link="logit"))
## boundary (singular) fit: see ?isSingular
summary(g24New)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: cbind(Te, Total) ~ Treatment + (1 | Block2)
##    Data: pair2_4_new
## 
##      AIC      BIC   logLik deviance df.resid 
##   3788.1   3829.2  -1886.1   3772.1     1242 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2080 -0.3468  0.0000  0.1375  5.9739 
## 
## Random effects:
##  Groups Name        Variance  Std.Dev. 
##  Block2 (Intercept) 8.736e-18 2.956e-09
## Number of obs: 1250, groups:  Block2, 2
## 
## Fixed effects:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -0.46012    0.03496 -13.163  < 2e-16 ***
## Treatment10Te10Tu48h -0.65741    0.05938 -11.071  < 2e-16 ***
## Treatment10Te48h10Tu  0.28270    0.04772   5.924 3.14e-09 ***
## Treatment19Te1Tu      0.43806    0.04483   9.771  < 2e-16 ***
## Treatment19Te48h1Tu   0.43257    0.04180  10.350  < 2e-16 ***
## Treatment1Te19Tu     -0.78099    0.08384  -9.316  < 2e-16 ***
## Treatment1Te19Tu48h  -2.52713    0.12192 -20.728  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) T10T10 T10T48 T19T1T T19T48 Tr1T19T
## Trt10T10T48 -0.589                                    
## Trt10T4810T -0.733  0.431                             
## Trtmnt19T1T -0.780  0.459  0.571                      
## Trtm19T481T -0.836  0.492  0.613  0.652               
## Trtmnt1T19T -0.417  0.245  0.305  0.325  0.349        
## Trtm1T19T48 -0.287  0.169  0.210  0.224  0.240  0.120 
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
Anova(g24New, type=3)
test_intOld<-testInteractions(g24Old, pairwise="Treatment", adjustment = "fdr")
test_intOld
test_intNew<-testInteractions(g24New, pairwise="Treatment", adjustment = "fdr")
test_intNew
#Now we just select the ones we want to compare
#effect of timing
test_intOld[c(1, 2,16, 21, 7 ),]
test_intNew[c(1, 2,16, 21, 7 ),]
#effect of density
test_intOld[c(3,5,11, 13,17 ),]
test_intNew[c(3,5,11, 13,17 ),]
test_intOld2<-testInteractions(g24Old, custom=list(Treatment=mat_overall), adjustment = "fdr")
test_intNew2<-testInteractions(g24New, custom=list(Treatment=mat_overall), adjustment = "fdr")

test_intOld2
test_intNew2

We have the same results as before! Which indicates that adding the different pairs of leaves first does not affect the impact of order of arrival or initial frequency.

3 - Impact of spatial variation

3.1 - Is spatial variation mediating coexistence

Analysing differences to the occupancy of single treatments

This analyses do not take into account leaf age!

Here we will calculate the percentage of individuals from each box that are in each leaf for the single regime. Then we will don the same for the treatments and will then compare the difference

Fig 3, S3

# Start with calculating total number of individuals per leaf
single_df<-subset(df, Treatment=="20Te" | Treatment=="20Tu")
str(single_df)
## 'data.frame':    714 obs. of  17 variables:
##  $ Box         : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Treatment   : chr  "20Te" "20Te" "20Te" "20Te" ...
##  $ Block       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Leaf        : int  2 2 2 2 2 2 3 3 3 3 ...
##  $ Leaflet     : int  1 2 3 1 2 3 1 2 3 4 ...
##  $ Direction   : Factor w/ 2 levels "Down","Up": 1 1 1 2 2 2 1 1 1 1 ...
##  $ Te          : int  2 20 0 13 10 0 7 15 2 4 ...
##  $ Tu          : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ leaf_age    : chr  "old" "old" "old" "old" ...
##  $ Density     : chr  "20:0" "20:0" "20:0" "20:0" ...
##  $ Timing      : chr  "single" "single" "single" "single" ...
##  $ ratio       : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Total       : int  2 20 0 13 10 0 7 15 2 4 ...
##  $ proportionTe: num  1 1 0 1 1 0 1 1 1 1 ...
##  $ proportionTu: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Block2      : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Box2        : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
total_box_single_leaf<-single_df %>%
  group_by(Treatment, Box, Leaf, Block) %>%
  summarize(sumTe=sum(Te, na.rm=TRUE),sumTu=sum(Tu, na.rm=TRUE))%>%
  as.data.frame()
## `summarise()` has grouped output by 'Treatment', 'Box', 'Leaf'. You can override using the `.groups` argument.
# Same for the treatments
total_box_treat_leaf<-df %>%
  group_by(Treatment, Box, Leaf, Block) %>%
  summarize(sumTe=sum(Te, na.rm=TRUE),sumTu=sum(Tu, na.rm=TRUE))%>%
  as.data.frame()
## `summarise()` has grouped output by 'Treatment', 'Box', 'Leaf'. You can override using the `.groups` argument.
## Calculating percentage

#For single
total_box_single_leaf$PercentTe<-sapply(c(1:length(total_box_single_leaf[,1])), function(x){
  
  aux<-subset(total_box_single_leaf, Treatment==total_box_single_leaf$Treatment[x] & Box==total_box_single_leaf$Box[x])
  
  aux_sum<-sum(aux$sumTe[1:4], na.rm=TRUE)
  
  a<-total_box_single_leaf$sumTe[x]/aux_sum
  
  if(a=="NaN"){
    a<-NA
  }
  
  a
})

total_box_single_leaf$PercentTu<-sapply(c(1:length(total_box_single_leaf[,1])), function(x){
  
  aux<-subset(total_box_single_leaf, Treatment==total_box_single_leaf$Treatment[x] & Box==total_box_single_leaf$Box[x])
  
  aux_sum<-sum(aux$sumTu[1:4], na.rm=TRUE)
  
  a<-total_box_single_leaf$sumTu[x]/aux_sum
  
  if(a=="NaN"){
    a<-NA
  }
  
  a
})

# For treatments
total_box_treat_leaf$PercentTe<-sapply(c(1:length(total_box_treat_leaf[,1])), function(x){
  
  aux<-subset(total_box_treat_leaf, Treatment==total_box_treat_leaf$Treatment[x] & Box==total_box_treat_leaf$Box[x])
  
  aux_sum<-sum(aux$sumTe[1:4], na.rm=TRUE)
  
  a<-total_box_treat_leaf$sumTe[x]/aux_sum
  
  if(a=="NaN"){
    a<-NA
  }
  
  a
})

total_box_treat_leaf$PercentTu<-sapply(c(1:length(total_box_treat_leaf[,1])), function(x){
  
  aux<-subset(total_box_treat_leaf, Treatment==total_box_treat_leaf$Treatment[x] & Box==total_box_treat_leaf$Box[x])
  
  aux_sum<-sum(aux$sumTu[1:4], na.rm=TRUE)
  
  a<-total_box_treat_leaf$sumTu[x]/aux_sum
  
  if(a=="NaN"){
    a<-NA
  }
  
  a
})

total_box_treat_leaf
### Now calculating the difference to the average percentage of the single regime

# first creating average data base for the single regimes

mean_percentage_single<-total_box_single_leaf %>%
  group_by(Treatment, Leaf)%>%
  summarize(meanTe=mean(PercentTe), meanTu=mean(PercentTu))%>%
  as.data.frame()
## `summarise()` has grouped output by 'Treatment'. You can override using the `.groups` argument.
mean_percentage_single
#Now adding to thr percent data frame
aux_diffTe<-t(sapply(c(1:length(total_box_treat_leaf[,1])), function(x){
  aux<-subset(mean_percentage_single, Leaf==total_box_treat_leaf$Leaf[x])
  
  te<-total_box_treat_leaf$PercentTe[x]-aux$meanTe[1]
  tu<-total_box_treat_leaf$PercentTu[x]-aux$meanTu[2] #because the first is always NA
  
  c(te, tu)
}))

colnames(aux_diffTe)<-c("diffTe", "diffTu")


total_box_treat_leaf<-as.data.frame(cbind(total_box_treat_leaf[,c(1:8)], aux_diffTe))
str(total_box_treat_leaf)
## 'data.frame':    360 obs. of  10 variables:
##  $ Treatment: chr  "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
##  $ Box      : int  1 1 1 1 2 2 2 2 3 3 ...
##  $ Leaf     : int  2 3 4 5 2 3 4 5 2 3 ...
##  $ Block    : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ sumTe    : int  22 22 80 8 68 39 102 28 19 22 ...
##  $ sumTu    : int  27 51 27 23 30 34 54 1 24 32 ...
##  $ PercentTe: num  0.1667 0.1667 0.6061 0.0606 0.2869 ...
##  $ PercentTu: num  0.211 0.398 0.211 0.18 0.252 ...
##  $ diffTe   : num  0.0459 -0.1074 0.155 -0.0935 0.1662 ...
##  $ diffTu   : num  0.0436 0.1188 -0.0968 -0.0656 0.0848 ...
total_box_treat_leaf$Treatment2<-factor(total_box_treat_leaf$Treatment, levels=c("1Te19Tu", "10Te10Tu","19Te1Tu", "1Te19Tu48h","10Te10Tu48h","10Te48h10Tu","19Te48h1Tu"))

total_box_treat_leaf$Treatment3<-mapvalues(total_box_treat_leaf$Treatment,c("1Te19Tu", "10Te10Tu","19Te1Tu", "1Te19Tu48h","10Te10Tu48h","10Te48h10Tu","19Te48h1Tu"), c("1:19\nsame time", "10:10\nsame time", "19:1\nsame time", "1:19\nTu first", "10:10\nTu first","10:10\nTe first", "19:1\nTe first") )

total_box_treat_leaf$Treatment4<-factor(total_box_treat_leaf$Treatment3, levels=c("1:19\nsame time", "10:10\nsame time", "19:1\nsame time", "1:19\nTu first", "10:10\nTu first","10:10\nTe first", "19:1\nTe first") )


plot_percent_Te<-ggplot(subset(total_box_treat_leaf,Treatment!="20Te" &  Treatment!="20Tu"), aes(x=as.factor(Leaf), y=diffTe, fill=Treatment4))+
  facet_grid(.~Treatment4)+
  geom_hline(yintercept = 0)+
  geom_boxplot()+
  theme_ines+
  scale_fill_manual(values = c("#e0f3f8", "#91bfdb", "#4575b4",  "#fee090", "#ffffbf","#d73027", "#fc8d59"), name="Treatment")+
  ylab(c("Percentage of T. evansi per leaf\nrelative to single control"))+
  xlab(c("Leaf Position"))+
  #scale_x_discrete(labels=c())+
  theme(legend.title = element_text(size=16, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
plot_percent_Te

save_plot("./Plots/Percent_Single_Te.pdf", width=25, height=15)

plot_percent_Tu<-ggplot(subset(total_box_treat_leaf,Treatment!="20Te" &  Treatment!="20Tu"), aes(x=as.factor(Leaf), y=diffTu, fill=Treatment4))+
  facet_grid(.~Treatment4)+
  geom_hline(yintercept = 0)+
  geom_boxplot()+
  theme_ines+
  scale_fill_manual(values = c("#e0f3f8", "#91bfdb", "#4575b4",  "#fee090", "#ffffbf","#d73027", "#fc8d59"), name="Treatment")+
  ylab(c("Density of T. urticae per leaf\nrelative to single control"))+
  xlab(c("Leaf Position"))+
  #scale_x_discrete(labels=c())+
  theme(legend.title = element_text(size=16, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
plot_percent_Tu
## Warning: Removed 4 rows containing non-finite values (stat_boxplot).

save_plot("./Plots/Percent_Single_Tu.pdf", width=25, height=15)
## Warning: Removed 4 rows containing non-finite values (stat_boxplot).
### For single

te_single_perc<-ggplot(subset(total_box_treat_leaf, Treatment=="20Te"), aes(x=as.factor(Leaf), y=PercentTe))+
  geom_boxplot(fill="#d73027")+
  theme_ines+
  ylab(c("Number Te females"))+
  xlab(c("Leaf Position"))+
  #scale_x_discrete(labels=c())+
  theme(legend.title = element_text(size=16, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
te_single_perc

save_plot("./Plots/Te_single_dist_perc.pdf", width=7.5, height=7.5)


tu_single_perc<-ggplot(subset(total_box_treat_leaf, Treatment=="20Tu"), aes(x=as.factor(Leaf), y=PercentTu))+
  geom_boxplot(fill="#fee090")+
  theme_ines+
  ylab(c("Number Tu females"))+
  xlab(c("Leaf Position"))+
  #scale_x_discrete(labels=c())+
  theme(legend.title = element_text(size=16, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
tu_single_perc

save_plot("./Plots/Tu_single_dist_perc.pdf", width=7.5, height=7.5)

Te prefers leaves 3 and 4 when alone and Tu does not show a strong preference for any leaf. When Tu arrives first, there is a decrease in Te occupancy for leaves 3 and 4, with a corresponding Tu occupancy of these leaves.

Stats percentage

total_box_treat_leaf$TreatmentTe<-factor(as.factor(total_box_treat_leaf$Treatment), levels=c("20Te","10Te10Tu", "10Te10Tu48h","10Te48h10Tu","19Te1Tu" ,"19Te48h1Tu", "1Te19Tu", "1Te19Tu48h","20Tu" ))
total_box_treat_leaf$TreatmentTu<-factor(as.factor(total_box_treat_leaf$Treatment), levels=c("20Tu","10Te10Tu", "10Te10Tu48h","10Te48h10Tu","19Te1Tu" ,"19Te48h1Tu", "1Te19Tu", "1Te19Tu48h","20Te" ))

total_box_treat_leaf$Block<-as.factor(total_box_treat_leaf$Block)
total_box_treat_leaf$Leaf2<-as.factor(total_box_treat_leaf$Leaf)

str(total_box_treat_leaf)
## 'data.frame':    360 obs. of  16 variables:
##  $ Treatment  : chr  "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
##  $ Box        : int  1 1 1 1 2 2 2 2 3 3 ...
##  $ Leaf       : int  2 3 4 5 2 3 4 5 2 3 ...
##  $ Block      : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ sumTe      : int  22 22 80 8 68 39 102 28 19 22 ...
##  $ sumTu      : int  27 51 27 23 30 34 54 1 24 32 ...
##  $ PercentTe  : num  0.1667 0.1667 0.6061 0.0606 0.2869 ...
##  $ PercentTu  : num  0.211 0.398 0.211 0.18 0.252 ...
##  $ diffTe     : num  0.0459 -0.1074 0.155 -0.0935 0.1662 ...
##  $ diffTu     : num  0.0436 0.1188 -0.0968 -0.0656 0.0848 ...
##  $ Treatment2 : Factor w/ 7 levels "1Te19Tu","10Te10Tu",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Treatment3 : chr  "10:10\nsame time" "10:10\nsame time" "10:10\nsame time" "10:10\nsame time" ...
##  $ Treatment4 : Factor w/ 7 levels "1:19\nsame time",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ TreatmentTe: Factor w/ 9 levels "20Te","10Te10Tu",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ TreatmentTu: Factor w/ 9 levels "20Tu","10Te10Tu",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Leaf2      : Factor w/ 4 levels "2","3","4","5": 1 2 3 4 1 2 3 4 1 2 ...
# adding the total mites
Total_mites<-t(sapply(c(1:length(total_box_treat_leaf[,1])), function(x){
  aux<-subset(total_box_treat_leaf, Treatment==total_box_treat_leaf$Treatment[x] & Box==total_box_treat_leaf$Box[x])
  
  te<-sum(aux$sumTe, na.rm = TRUE)
  tu<-sum(aux$sumTu, na.rm = TRUE)
  
  c(te, tu)
}))

colnames(Total_mites)<-c("TotalTe", "TotalTu")

total_box_treat_leaf<-as.data.frame(cbind(total_box_treat_leaf[, c(1:16)], Total_mites))


forTe<-subset(total_box_treat_leaf, Treatment!="20Tu")
forTe<-droplevels(forTe)

forTu<-subset(total_box_treat_leaf, Treatment!="20Te")
forTu<-droplevels(forTu)

str(forTe)
## 'data.frame':    320 obs. of  18 variables:
##  $ Treatment  : chr  "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
##  $ Box        : int  1 1 1 1 2 2 2 2 3 3 ...
##  $ Leaf       : int  2 3 4 5 2 3 4 5 2 3 ...
##  $ Block      : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ sumTe      : int  22 22 80 8 68 39 102 28 19 22 ...
##  $ sumTu      : int  27 51 27 23 30 34 54 1 24 32 ...
##  $ PercentTe  : num  0.1667 0.1667 0.6061 0.0606 0.2869 ...
##  $ PercentTu  : num  0.211 0.398 0.211 0.18 0.252 ...
##  $ diffTe     : num  0.0459 -0.1074 0.155 -0.0935 0.1662 ...
##  $ diffTu     : num  0.0436 0.1188 -0.0968 -0.0656 0.0848 ...
##  $ Treatment2 : Factor w/ 7 levels "1Te19Tu","10Te10Tu",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Treatment3 : chr  "10:10\nsame time" "10:10\nsame time" "10:10\nsame time" "10:10\nsame time" ...
##  $ Treatment4 : Factor w/ 7 levels "1:19\nsame time",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ TreatmentTe: Factor w/ 8 levels "20Te","10Te10Tu",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ TreatmentTu: Factor w/ 8 levels "10Te10Tu","10Te10Tu48h",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Leaf2      : Factor w/ 4 levels "2","3","4","5": 1 2 3 4 1 2 3 4 1 2 ...
##  $ TotalTe    : int  132 132 132 132 237 237 237 237 117 117 ...
##  $ TotalTu    : int  128 128 128 128 119 119 119 119 133 133 ...
## using cbind

# when performing the glmer the block shows variance 0, so I will run the model without the random effect
perc_te<-glm(cbind(sumTe, TotalTe)~TreatmentTe*Leaf2, data=forTe, family=binomial) 
summary(perc_te)
## 
## Call:
## glm(formula = cbind(sumTe, TotalTe) ~ TreatmentTe * Leaf2, family = binomial, 
##     data = forTe)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -14.8605   -3.1028   -0.6477    1.7965   19.5147  
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   -2.126712   0.050493 -42.119  < 2e-16 ***
## TreatmentTe10Te10Tu            0.264994   0.075427   3.513 0.000443 ***
## TreatmentTe10Te10Tu48h         0.195713   0.087790   2.229 0.025793 *  
## TreatmentTe10Te48h10Tu         0.175420   0.069427   2.527 0.011514 *  
## TreatmentTe19Te1Tu             0.003126   0.064908   0.048 0.961594    
## TreatmentTe19Te48h1Tu          0.292952   0.060735   4.823 1.41e-06 ***
## TreatmentTe1Te19Tu            -0.121945   0.173747  -0.702 0.482771    
## TreatmentTe1Te19Tu48h          0.845778   0.175911   4.808 1.52e-06 ***
## Leaf23                         0.882468   0.061348  14.385  < 2e-16 ***
## Leaf24                         1.304447   0.058644  22.243  < 2e-16 ***
## Leaf25                         0.250555   0.067780   3.697 0.000219 ***
## TreatmentTe10Te10Tu:Leaf23    -0.629343   0.097135  -6.479 9.23e-11 ***
## TreatmentTe10Te10Tu48h:Leaf23 -0.096881   0.107832  -0.898 0.368951    
## TreatmentTe10Te48h10Tu:Leaf23 -0.198307   0.085555  -2.318 0.020455 *  
## TreatmentTe19Te1Tu:Leaf23     -0.092391   0.079242  -1.166 0.243639    
## TreatmentTe19Te48h1Tu:Leaf23  -0.245901   0.074697  -3.292 0.000995 ***
## TreatmentTe1Te19Tu:Leaf23     -0.308667   0.219374  -1.407 0.159417    
## TreatmentTe1Te19Tu48h:Leaf23  -1.051544   0.254258  -4.136 3.54e-05 ***
## TreatmentTe10Te10Tu:Leaf24    -0.176381   0.088773  -1.987 0.046936 *  
## TreatmentTe10Te10Tu48h:Leaf24 -0.591232   0.107040  -5.523 3.32e-08 ***
## TreatmentTe10Te48h10Tu:Leaf24 -0.293750   0.081950  -3.585 0.000338 ***
## TreatmentTe19Te1Tu:Leaf24     -0.135555   0.075777  -1.789 0.073636 .  
## TreatmentTe19Te48h1Tu:Leaf24  -0.453507   0.071798  -6.316 2.68e-10 ***
## TreatmentTe1Te19Tu:Leaf24      0.269364   0.197225   1.366 0.172010    
## TreatmentTe1Te19Tu48h:Leaf24  -1.239908   0.242570  -5.112 3.20e-07 ***
## TreatmentTe10Te10Tu:Leaf25    -0.195058   0.103570  -1.883 0.059654 .  
## TreatmentTe10Te10Tu48h:Leaf25  0.257564   0.114521   2.249 0.024508 *  
## TreatmentTe10Te48h10Tu:Leaf25  0.018060   0.093096   0.194 0.846182    
## TreatmentTe19Te1Tu:Leaf25      0.411658   0.084863   4.851 1.23e-06 ***
## TreatmentTe19Te48h1Tu:Leaf25  -0.225335   0.082754  -2.723 0.006470 ** 
## TreatmentTe1Te19Tu:Leaf25      0.378054   0.219555   1.722 0.085086 .  
## TreatmentTe1Te19Tu48h:Leaf25  -0.623230   0.267245  -2.332 0.019698 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 9140.8  on 319  degrees of freedom
## Residual deviance: 5961.6  on 288  degrees of freedom
## AIC: 7590.8
## 
## Number of Fisher Scoring iterations: 5
perc_tu<-glm(cbind(sumTu, TotalTu)~TreatmentTu*Leaf2, data=forTu, family=binomial) 

summary(perc_te)
## 
## Call:
## glm(formula = cbind(sumTe, TotalTe) ~ TreatmentTe * Leaf2, family = binomial, 
##     data = forTe)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -14.8605   -3.1028   -0.6477    1.7965   19.5147  
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   -2.126712   0.050493 -42.119  < 2e-16 ***
## TreatmentTe10Te10Tu            0.264994   0.075427   3.513 0.000443 ***
## TreatmentTe10Te10Tu48h         0.195713   0.087790   2.229 0.025793 *  
## TreatmentTe10Te48h10Tu         0.175420   0.069427   2.527 0.011514 *  
## TreatmentTe19Te1Tu             0.003126   0.064908   0.048 0.961594    
## TreatmentTe19Te48h1Tu          0.292952   0.060735   4.823 1.41e-06 ***
## TreatmentTe1Te19Tu            -0.121945   0.173747  -0.702 0.482771    
## TreatmentTe1Te19Tu48h          0.845778   0.175911   4.808 1.52e-06 ***
## Leaf23                         0.882468   0.061348  14.385  < 2e-16 ***
## Leaf24                         1.304447   0.058644  22.243  < 2e-16 ***
## Leaf25                         0.250555   0.067780   3.697 0.000219 ***
## TreatmentTe10Te10Tu:Leaf23    -0.629343   0.097135  -6.479 9.23e-11 ***
## TreatmentTe10Te10Tu48h:Leaf23 -0.096881   0.107832  -0.898 0.368951    
## TreatmentTe10Te48h10Tu:Leaf23 -0.198307   0.085555  -2.318 0.020455 *  
## TreatmentTe19Te1Tu:Leaf23     -0.092391   0.079242  -1.166 0.243639    
## TreatmentTe19Te48h1Tu:Leaf23  -0.245901   0.074697  -3.292 0.000995 ***
## TreatmentTe1Te19Tu:Leaf23     -0.308667   0.219374  -1.407 0.159417    
## TreatmentTe1Te19Tu48h:Leaf23  -1.051544   0.254258  -4.136 3.54e-05 ***
## TreatmentTe10Te10Tu:Leaf24    -0.176381   0.088773  -1.987 0.046936 *  
## TreatmentTe10Te10Tu48h:Leaf24 -0.591232   0.107040  -5.523 3.32e-08 ***
## TreatmentTe10Te48h10Tu:Leaf24 -0.293750   0.081950  -3.585 0.000338 ***
## TreatmentTe19Te1Tu:Leaf24     -0.135555   0.075777  -1.789 0.073636 .  
## TreatmentTe19Te48h1Tu:Leaf24  -0.453507   0.071798  -6.316 2.68e-10 ***
## TreatmentTe1Te19Tu:Leaf24      0.269364   0.197225   1.366 0.172010    
## TreatmentTe1Te19Tu48h:Leaf24  -1.239908   0.242570  -5.112 3.20e-07 ***
## TreatmentTe10Te10Tu:Leaf25    -0.195058   0.103570  -1.883 0.059654 .  
## TreatmentTe10Te10Tu48h:Leaf25  0.257564   0.114521   2.249 0.024508 *  
## TreatmentTe10Te48h10Tu:Leaf25  0.018060   0.093096   0.194 0.846182    
## TreatmentTe19Te1Tu:Leaf25      0.411658   0.084863   4.851 1.23e-06 ***
## TreatmentTe19Te48h1Tu:Leaf25  -0.225335   0.082754  -2.723 0.006470 ** 
## TreatmentTe1Te19Tu:Leaf25      0.378054   0.219555   1.722 0.085086 .  
## TreatmentTe1Te19Tu48h:Leaf25  -0.623230   0.267245  -2.332 0.019698 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 9140.8  on 319  degrees of freedom
## Residual deviance: 5961.6  on 288  degrees of freedom
## AIC: 7590.8
## 
## Number of Fisher Scoring iterations: 5
Anova(perc_te)
summary(perc_tu)
## 
## Call:
## glm(formula = cbind(sumTu, TotalTu) ~ TreatmentTu * Leaf2, family = binomial, 
##     data = forTu)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -16.6458   -2.2388   -0.5321    1.7716   14.2185  
## 
## Coefficients:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   -1.79096    0.07472 -23.970  < 2e-16 ***
## TreatmentTu10Te10Tu           -0.09638    0.10333  -0.933 0.350990    
## TreatmentTu10Te10Tu48h        -0.32784    0.10068  -3.256 0.001129 ** 
## TreatmentTu10Te48h10Tu         0.05277    0.12066   0.437 0.661830    
## TreatmentTu19Te1Tu            -0.11858    0.27810  -0.426 0.669824    
## TreatmentTu19Te48h1Tu         -1.67997    0.42122  -3.988 6.65e-05 ***
## TreatmentTu1Te19Tu             0.26059    0.10530   2.475 0.013330 *  
## TreatmentTu1Te19Tu48h         -0.13117    0.08980  -1.461 0.144105    
## Leaf23                         0.57117    0.09531   5.993 2.06e-09 ***
## Leaf24                         0.54654    0.09566   5.713 1.11e-08 ***
## Leaf25                         0.40387    0.09786   4.127 3.67e-05 ***
## TreatmentTu10Te10Tu:Leaf23     0.06039    0.13116   0.460 0.645221    
## TreatmentTu10Te10Tu48h:Leaf23 -0.02677    0.12839  -0.208 0.834854    
## TreatmentTu10Te48h10Tu:Leaf23  0.16676    0.15183   1.098 0.272043    
## TreatmentTu19Te1Tu:Leaf23      0.21159    0.34449   0.614 0.539074    
## TreatmentTu19Te48h1Tu:Leaf23   2.22204    0.44307   5.015 5.30e-07 ***
## TreatmentTu1Te19Tu:Leaf23     -0.43990    0.13977  -3.147 0.001648 ** 
## TreatmentTu1Te19Tu48h:Leaf23   0.16020    0.11368   1.409 0.158770    
## TreatmentTu10Te10Tu:Leaf24    -0.05358    0.13274  -0.404 0.686488    
## TreatmentTu10Te10Tu48h:Leaf24  0.78653    0.12355   6.366 1.94e-10 ***
## TreatmentTu10Te48h10Tu:Leaf24  0.03358    0.15409   0.218 0.827483    
## TreatmentTu19Te1Tu:Leaf24      0.39444    0.33847   1.165 0.243868    
## TreatmentTu19Te48h1Tu:Leaf24   1.57372    0.45407   3.466 0.000529 ***
## TreatmentTu1Te19Tu:Leaf24     -0.14865    0.13665  -1.088 0.276694    
## TreatmentTu1Te19Tu48h:Leaf24   0.26005    0.11364   2.288 0.022122 *  
## TreatmentTu10Te10Tu:Leaf25     0.33045    0.13215   2.501 0.012400 *  
## TreatmentTu10Te10Tu48h:Leaf25  0.18662    0.12988   1.437 0.150763    
## TreatmentTu10Te48h10Tu:Leaf25 -0.61563    0.17125  -3.595 0.000325 ***
## TreatmentTu19Te1Tu:Leaf25     -0.40387    0.39127  -1.032 0.301981    
## TreatmentTu19Te48h1Tu:Leaf25   1.46793    0.46070   3.186 0.001441 ** 
## TreatmentTu1Te19Tu:Leaf25     -0.41296    0.14362  -2.875 0.004035 ** 
## TreatmentTu1Te19Tu48h:Leaf25   0.01307    0.11747   0.111 0.911384    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 5399.3  on 315  degrees of freedom
## Residual deviance: 4469.2  on 284  degrees of freedom
## AIC: 5776.8
## 
## Number of Fisher Scoring iterations: 5
Anova(perc_tu)
testInteractions(perc_te, fixed = c("Leaf2"), pairwise = c("TreatmentTe"),adjustment = "fdr")
testInteractions(perc_tu, fixed = c("Leaf2"), pairwise = c("TreatmentTu"),adjustment = "fdr")

Fig S4

total_box_single_24old<-subset(single_df, (Leaf==2 & leaf_age=="old") | (Leaf==4 & leaf_age=="old") | (Leaf==3 & leaf_age=="new") | (Leaf==5 & leaf_age=="new")) %>%
  group_by(Treatment, Box, Leaf, Block) %>%
  summarize(sumTe=sum(Te, na.rm=TRUE),sumTu=sum(Tu, na.rm=TRUE))%>%
  as.data.frame()
## `summarise()` has grouped output by 'Treatment', 'Box', 'Leaf'. You can override using the `.groups` argument.
total_box_single_24new<-subset(single_df, (Leaf==3 & leaf_age=="old") | (Leaf==5 & leaf_age=="old") | (Leaf==2 & leaf_age=="new") | (Leaf==4 & leaf_age=="new")) %>%
  group_by(Treatment, Box, Leaf, Block) %>%
  summarize(sumTe=sum(Te, na.rm=TRUE),sumTu=sum(Tu, na.rm=TRUE))%>%
  as.data.frame()
## `summarise()` has grouped output by 'Treatment', 'Box', 'Leaf'. You can override using the `.groups` argument.
# Same for the treatments
total_box_treat_24old<-subset(df_wsingle, (Leaf==2 & leaf_age=="old") | (Leaf==4 & leaf_age=="old") | (Leaf==3 & leaf_age=="new") | (Leaf==5 & leaf_age=="new")) %>%
  group_by(Treatment, Box, Leaf, Block) %>%
  summarize(sumTe=sum(Te, na.rm=TRUE),sumTu=sum(Tu, na.rm=TRUE))%>%
  as.data.frame()
## `summarise()` has grouped output by 'Treatment', 'Box', 'Leaf'. You can override using the `.groups` argument.
total_box_treat_24new<-subset(df_wsingle, (Leaf==3 & leaf_age=="old") | (Leaf==5 & leaf_age=="old") | (Leaf==2 & leaf_age=="new") | (Leaf==4 & leaf_age=="new")) %>%
  group_by(Treatment, Box, Leaf, Block) %>%
  summarize(sumTe=sum(Te, na.rm=TRUE),sumTu=sum(Tu, na.rm=TRUE))%>%
  as.data.frame()
## `summarise()` has grouped output by 'Treatment', 'Box', 'Leaf'. You can override using the `.groups` argument.
## Calculating percentage

#For single
total_box_single_24old$PercentTe<-sapply(c(1:length(total_box_single_24old[,1])), function(x){
  
  aux<-subset(total_box_single_24old, Treatment==total_box_single_24old$Treatment[x] & Box==total_box_single_24old$Box[x])
  
  aux_sum<-sum(aux$sumTe[1:4], na.rm=TRUE)
  
  a<-total_box_single_24old$sumTe[x]/aux_sum
  
  if(a=="NaN"){
    a<-NA
  }
  
  a
})

total_box_single_24old$PercentTu<-sapply(c(1:length(total_box_single_24old[,1])), function(x){
  
  aux<-subset(total_box_single_24old, Treatment==total_box_single_24old$Treatment[x] & Box==total_box_single_24old$Box[x])
  
  aux_sum<-sum(aux$sumTu[1:4], na.rm=TRUE)
  
  a<-total_box_single_24old$sumTu[x]/aux_sum
  
  if(a=="NaN"){
    a<-NA
  }
  
  a
})


total_box_single_24new$PercentTe<-sapply(c(1:length(total_box_single_24new[,1])), function(x){
  
  aux<-subset(total_box_single_24new, Treatment==total_box_single_24new$Treatment[x] & Box==total_box_single_24new$Box[x])
  
  aux_sum<-sum(aux$sumTe[1:4], na.rm=TRUE)
  
  a<-total_box_single_24new$sumTe[x]/aux_sum
  
  if(a=="NaN"){
    a<-NA
  }
  
  a
})

total_box_single_24new$PercentTu<-sapply(c(1:length(total_box_single_24new[,1])), function(x){
  
  aux<-subset(total_box_single_24new, Treatment==total_box_single_24new$Treatment[x] & Box==total_box_single_24new$Box[x])
  
  aux_sum<-sum(aux$sumTu[1:4], na.rm=TRUE)
  
  a<-total_box_single_24new$sumTu[x]/aux_sum
  
  if(a=="NaN"){
    a<-NA
  }
  
  a
})

# For treatments
total_box_treat_24old$PercentTe<-sapply(c(1:length(total_box_treat_24old[,1])), function(x){
  
  aux<-subset(total_box_treat_24old, Treatment==total_box_treat_24old$Treatment[x] & Box==total_box_treat_24old$Box[x])
  
  aux_sum<-sum(aux$sumTe[1:4], na.rm=TRUE)
  
  a<-total_box_treat_24old$sumTe[x]/aux_sum
  
  if(a=="NaN"){
    a<-NA
  }
  
  a
})

total_box_treat_24old$PercentTu<-sapply(c(1:length(total_box_treat_24old[,1])), function(x){
  
  aux<-subset(total_box_treat_24old, Treatment==total_box_treat_24old$Treatment[x] & Box==total_box_treat_24old$Box[x])
  
  aux_sum<-sum(aux$sumTu[1:4], na.rm=TRUE)
  
  a<-total_box_treat_24old$sumTu[x]/aux_sum
  
  if(a=="NaN"){
    a<-NA
  }
  
  a
})

total_box_treat_24new$PercentTe<-sapply(c(1:length(total_box_treat_24new[,1])), function(x){
  
  aux<-subset(total_box_treat_24new, Treatment==total_box_treat_24new$Treatment[x] & Box==total_box_treat_24new$Box[x])
  
  aux_sum<-sum(aux$sumTe[1:4], na.rm=TRUE)
  
  a<-total_box_treat_24new$sumTe[x]/aux_sum
  
  if(a=="NaN"){
    a<-NA
  }
  
  a
})

total_box_treat_24new$PercentTu<-sapply(c(1:length(total_box_treat_24new[,1])), function(x){
  
  aux<-subset(total_box_treat_24new, Treatment==total_box_treat_24new$Treatment[x] & Box==total_box_treat_24new$Box[x])
  
  aux_sum<-sum(aux$sumTu[1:4], na.rm=TRUE)
  
  a<-total_box_treat_24new$sumTu[x]/aux_sum
  
  if(a=="NaN"){
    a<-NA
  }
  
  a
})


### Now calculating the difference to the average percentage of the single regime

# first creating average data base for the single regimes

mean_percentage_single_24old<-total_box_single_24old %>%
  group_by(Treatment, Leaf)%>%
  summarize(meanTe=mean(PercentTe), meanTu=mean(PercentTu))%>%
  as.data.frame()
## `summarise()` has grouped output by 'Treatment'. You can override using the `.groups` argument.
mean_percentage_single_24new<-total_box_single_24new %>%
  group_by(Treatment, Leaf)%>%
  summarize(meanTe=mean(PercentTe), meanTu=mean(PercentTu))%>%
  as.data.frame()
## `summarise()` has grouped output by 'Treatment'. You can override using the `.groups` argument.
#Now adding to the percent data frame
aux_diffTe_24old<-t(sapply(c(1:length(total_box_treat_24old[,1])), function(x){
  aux<-subset(mean_percentage_single_24old, Leaf==total_box_treat_24old$Leaf[x])
  
  te<-total_box_treat_24old$PercentTe[x]-aux$meanTe[1]
  tu<-total_box_treat_24old$PercentTu[x]-aux$meanTu[2] #because the first is always NA
  
  c(te, tu)
}))

aux_diffTe_24new<-t(sapply(c(1:length(total_box_treat_24new[,1])), function(x){
  aux<-subset(mean_percentage_single_24new, Leaf==total_box_treat_24new$Leaf[x])
  
  te<-total_box_treat_24new$PercentTe[x]-aux$meanTe[1]
  tu<-total_box_treat_24new$PercentTu[x]-aux$meanTu[2] #because the first is always NA
  
  c(te, tu)
}))
colnames(aux_diffTe_24old)<-c("diffTe", "diffTu")

colnames(aux_diffTe_24new)<-c("diffTe", "diffTu")


total_box_treat_24old<-as.data.frame(cbind(total_box_treat_24old[,c(1:8)], aux_diffTe_24old))
str(total_box_treat_24old)
## 'data.frame':    140 obs. of  10 variables:
##  $ Treatment: chr  "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
##  $ Box      : int  1 1 1 1 3 3 3 3 5 5 ...
##  $ Leaf     : int  2 3 4 5 2 3 4 5 2 3 ...
##  $ Block    : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ sumTe    : int  22 22 80 8 19 22 44 32 77 54 ...
##  $ sumTu    : int  27 51 27 23 24 32 33 44 6 21 ...
##  $ PercentTe: num  0.1667 0.1667 0.6061 0.0606 0.1624 ...
##  $ PercentTu: num  0.211 0.398 0.211 0.18 0.18 ...
##  $ diffTe   : num  0.0264 -0.0446 0.109 -0.0908 0.0221 ...
##  $ diffTu   : num  0.14013 -0.00689 0.05903 -0.19226 0.10964 ...
total_box_treat_24new<-as.data.frame(cbind(total_box_treat_24new[,c(1:8)], aux_diffTe_24new))
str(total_box_treat_24new)
## 'data.frame':    140 obs. of  10 variables:
##  $ Treatment: chr  "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
##  $ Box      : int  2 2 2 2 4 4 4 4 6 6 ...
##  $ Leaf     : int  2 3 4 5 2 3 4 5 2 3 ...
##  $ Block    : int  1 1 1 1 1 1 1 1 2 2 ...
##  $ sumTe    : int  68 39 102 28 49 38 85 142 32 28 ...
##  $ sumTu    : int  30 34 54 1 34 18 45 26 20 13 ...
##  $ PercentTe: num  0.287 0.165 0.43 0.118 0.156 ...
##  $ PercentTu: num  0.2521 0.2857 0.4538 0.0084 0.2764 ...
##  $ diffTe   : num  0.1857 -0.1724 0.0254 -0.0387 0.0548 ...
##  $ diffTu   : num  -0.0117 0.1318 -0.0098 -0.1103 0.0126 ...
total_box_treat_24old$Treatment2<-factor(total_box_treat_24old$Treatment, levels=c("1Te19Tu", "10Te10Tu","19Te1Tu", "1Te19Tu48h","10Te10Tu48h","10Te48h10Tu","19Te48h1Tu"))

total_box_treat_24old$Treatment3<-mapvalues(total_box_treat_24old$Treatment,c("1Te19Tu", "10Te10Tu","19Te1Tu", "1Te19Tu48h","10Te10Tu48h","10Te48h10Tu","19Te48h1Tu"), c("1:19\nsame time", "10:10\nsame time", "19:1\nsame time", "1:19\nTu first", "10:10\nTu first","10:10\nTe first", "19:1\nTe first") )

total_box_treat_24old$Treatment4<-factor(total_box_treat_24old$Treatment3, levels=c("1:19\nsame time", "10:10\nsame time", "19:1\nsame time", "1:19\nTu first", "10:10\nTu first","10:10\nTe first", "19:1\nTe first") )

total_box_treat_24new$Treatment2<-factor(total_box_treat_24new$Treatment, levels=c("1Te19Tu", "10Te10Tu","19Te1Tu", "1Te19Tu48h","10Te10Tu48h","10Te48h10Tu","19Te48h1Tu"))

total_box_treat_24new$Treatment3<-mapvalues(total_box_treat_24new$Treatment,c("1Te19Tu", "10Te10Tu","19Te1Tu", "1Te19Tu48h","10Te10Tu48h","10Te48h10Tu","19Te48h1Tu"), c("1:19\nsame time", "10:10\nsame time", "19:1\nsame time", "1:19\nTu first", "10:10\nTu first","10:10\nTe first", "19:1\nTe first") )

total_box_treat_24new$Treatment4<-factor(total_box_treat_24new$Treatment3, levels=c("1:19\nsame time", "10:10\nsame time", "19:1\nsame time", "1:19\nTu first", "10:10\nTu first","10:10\nTe first", "19:1\nTe first") )

### For single

te_single_24old<-ggplot(subset(total_box_single_24old, Treatment=="20Te"), aes(x=as.factor(Leaf), y=PercentTe))+
  geom_boxplot(fill="#d73027")+
  theme_ines+
  ylab(c("Number Te females"))+
  xlab(c("Leaf Position"))+
  #scale_x_discrete(labels=c())+
  theme(legend.title = element_text(size=16, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
te_single_24old

save_plot("./Plots/Te_single_dist_24old.pdf", width=7.5, height=7.5)

te_single_24new<-ggplot(subset(total_box_single_24new, Treatment=="20Te"), aes(x=as.factor(Leaf), y=PercentTe))+
  geom_boxplot(fill="#d73027")+
  theme_ines+
  ylab(c("Number Te females"))+
  xlab(c("Leaf Position"))+
  #scale_x_discrete(labels=c())+
  theme(legend.title = element_text(size=16, face="bnew"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
te_single_24new

save_plot("./Plots/Te_single_dist_24new.pdf", width=7.5, height=7.5)


tu_single_24old<-ggplot(subset(total_box_single_24old, Treatment=="20Tu"), aes(x=as.factor(Leaf), y=PercentTu))+
  geom_boxplot(fill="#fee090")+
  theme_ines+
  ylab(c("Number Tu females"))+
  xlab(c("Leaf Position"))+
  ylim(c(0,1))+
  #scale_x_discrete(labels=c())+
  theme(legend.title = element_text(size=16, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
tu_single_24old

save_plot("./Plots/Tu_single_dist_24old.pdf", width=7.5, height=7.5)

tu_single_24new<-ggplot(subset(total_box_single_24new, Treatment=="20Tu"), aes(x=as.factor(Leaf), y=PercentTu))+
  geom_boxplot(fill="#fee090")+
  theme_ines+
  ylab(c("Number Tu females"))+
  xlab(c("Leaf Position"))+
  #scale_x_discrete(labels=c())+
  ylim(c(0,1))+
  theme(legend.title = element_text(size=16, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
tu_single_24new

save_plot("./Plots/Tu_single_dist_24new.pdf", width=7.5, height=7.5)



#### Single graph with both old and new BUT WITH THE DIFFERENCES TO THE RESPECTIVE CONTROLS

new_total_box<-as.data.frame(rbind(total_box_treat_24new, total_box_treat_24new))


one_diffplot_Te<-ggplot(subset(new_total_box,Treatment!="20Te" &  Treatment!="20Tu"), aes(x=as.factor(Leaf), y=diffTe, fill=Treatment4))+
  facet_grid(.~Treatment4)+
  geom_hline(yintercept = 0)+
  geom_boxplot()+
  theme_ines+
  scale_fill_manual(values = c("#e0f3f8", "#91bfdb", "#4575b4",  "#fee090", "#ffffbf","#d73027", "#fc8d59"), name="Treatment")+
  ylab(c("Percentage of T. evansi per leaf\nrelative to single control"))+
  xlab(c("Leaf Position"))+
  #scale_x_discrete(labels=c())+
  theme(legend.title = element_text(size=16, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
one_diffplot_Te

save_plot("./Plots/one_diffplot_Te.pdf", width=25, height=15)

one_diffplot_Tu<-ggplot(subset(new_total_box,Treatment!="20Tu" &  Treatment!="20Tu"), aes(x=as.factor(Leaf), y=diffTu, fill=Treatment4))+
  facet_grid(.~Treatment4)+
  geom_hline(yintercept = 0)+
  geom_boxplot()+
  theme_ines+
  scale_fill_manual(values = c("#e0f3f8", "#91bfdb", "#4575b4",  "#fee090", "#ffffbf","#d73027", "#fc8d59"), name="Treatment")+
  ylab(c("Percentage of T. urticae per leaf\nrelative to single control"))+
  xlab(c("Leaf Position"))+
  #scale_x_discreTu(labels=c())+
  theme(legend.title = element_text(size=16, face="bold"), axis.text.x = element_text(angle=20, vjust = 0.75), legend.position = "bottom")
one_diffplot_Tu

save_plot("./Plots/one_diffplot_Tu.pdf", width=25, height=15)

Statistics (accounting for leaf order)

Distribution

hist(new_total_box$diffTe)

hist(new_total_box$diffTu)

descdist(new_total_box$diffTe, discrete = FALSE, boot=500)

## summary statistics
## ------
## min:  -0.4049489   max:  0.5098335 
## median:  -0.01485167 
## mean:  6.147063e-18 
## estimated sd:  0.1639964 
## estimated skewness:  0.4349104 
## estimated kurtosis:  3.35588
descdist(new_total_box$diffTu, discrete = FALSE, boot=500)

## summary statistics
## ------
## min:  -0.463581   max:  0.8461161 
## median:  -0.03990551 
## mean:  -6.938894e-18 
## estimated sd:  0.2168845 
## estimated skewness:  1.133723 
## estimated kurtosis:  5.350042
#Log normal fits both

Statistics

str(new_total_box)
## 'data.frame':    280 obs. of  13 variables:
##  $ Treatment : chr  "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
##  $ Box       : int  2 2 2 2 4 4 4 4 6 6 ...
##  $ Leaf      : int  2 3 4 5 2 3 4 5 2 3 ...
##  $ Block     : int  1 1 1 1 1 1 1 1 2 2 ...
##  $ sumTe     : int  68 39 102 28 49 38 85 142 32 28 ...
##  $ sumTu     : int  30 34 54 1 34 18 45 26 20 13 ...
##  $ PercentTe : num  0.287 0.165 0.43 0.118 0.156 ...
##  $ PercentTu : num  0.2521 0.2857 0.4538 0.0084 0.2764 ...
##  $ diffTe    : num  0.1857 -0.1724 0.0254 -0.0387 0.0548 ...
##  $ diffTu    : num  -0.0117 0.1318 -0.0098 -0.1103 0.0126 ...
##  $ Treatment2: Factor w/ 7 levels "1Te19Tu","10Te10Tu",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Treatment3: chr  "10:10\nsame time" "10:10\nsame time" "10:10\nsame time" "10:10\nsame time" ...
##  $ Treatment4: Factor w/ 7 levels "1:19\nsame time",..: 2 2 2 2 2 2 2 2 2 2 ...
# Making leaf and Block as factors
new_total_box$Leaf2<-as.factor(new_total_box$Leaf)
new_total_box$Block2<-as.factor(new_total_box$Block)

mla_Te<-lmer(log(diffTe+1)~Treatment2*Leaf2+ (1|Block2), data=new_total_box)
## boundary (singular) fit: see ?isSingular
mla_Tu<-lmer(log(diffTu+1)~Treatment2*Leaf2+ (1|Block2), data=new_total_box)
## boundary (singular) fit: see ?isSingular
summary(mla_Te)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(diffTe + 1) ~ Treatment2 * Leaf2 + (1 | Block2)
##    Data: new_total_box
## 
## REML criterion at convergence: -172.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7410 -0.5897 -0.1344  0.5338  3.0181 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Block2   (Intercept) 0.00000  0.0000  
##  Residual             0.02289  0.1513  
## Number of obs: 280, groups:  Block2, 2
## 
## Fixed effects:
##                               Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)                   -0.04177    0.04784 252.00000  -0.873  0.38343   
## Treatment210Te10Tu             0.11006    0.06766 252.00000   1.627  0.10506   
## Treatment219Te1Tu              0.08407    0.06766 252.00000   1.243  0.21518   
## Treatment21Te19Tu48h           0.12896    0.06766 252.00000   1.906  0.05780 . 
## Treatment210Te10Tu48h          0.12091    0.06766 252.00000   1.787  0.07514 . 
## Treatment210Te48h10Tu          0.13183    0.06766 252.00000   1.948  0.05248 . 
## Treatment219Te48h1Tu           0.16474    0.06766 252.00000   2.435  0.01560 * 
## Leaf23                        -0.06911    0.06766 252.00000  -1.021  0.30800   
## Leaf24                         0.08866    0.06766 252.00000   1.310  0.19124   
## Leaf25                         0.05220    0.06766 252.00000   0.772  0.44111   
## Treatment210Te10Tu:Leaf23     -0.17565    0.09569 252.00000  -1.836  0.06759 . 
## Treatment219Te1Tu:Leaf23      -0.00421    0.09569 252.00000  -0.044  0.96494   
## Treatment21Te19Tu48h:Leaf23   -0.13432    0.09569 252.00000  -1.404  0.16164   
## Treatment210Te10Tu48h:Leaf23   0.03848    0.09569 252.00000   0.402  0.68791   
## Treatment210Te48h10Tu:Leaf23  -0.06526    0.09569 252.00000  -0.682  0.49586   
## Treatment219Te48h1Tu:Leaf23   -0.08879    0.09569 252.00000  -0.928  0.35436   
## Treatment210Te10Tu:Leaf24     -0.10078    0.09569 252.00000  -1.053  0.29326   
## Treatment219Te1Tu:Leaf24      -0.05353    0.09569 252.00000  -0.559  0.57635   
## Treatment21Te19Tu48h:Leaf24   -0.28975    0.09569 252.00000  -3.028  0.00272 **
## Treatment210Te10Tu48h:Leaf24  -0.25875    0.09569 252.00000  -2.704  0.00731 **
## Treatment210Te48h10Tu:Leaf24  -0.27732    0.09569 252.00000  -2.898  0.00408 **
## Treatment219Te48h1Tu:Leaf24   -0.27320    0.09569 252.00000  -2.855  0.00466 **
## Treatment210Te10Tu:Leaf25     -0.11585    0.09569 252.00000  -1.211  0.22716   
## Treatment219Te1Tu:Leaf25      -0.21379    0.09569 252.00000  -2.234  0.02634 * 
## Treatment21Te19Tu48h:Leaf25   -0.09856    0.09569 252.00000  -1.030  0.30396   
## Treatment210Te10Tu48h:Leaf25  -0.20139    0.09569 252.00000  -2.105  0.03631 * 
## Treatment210Te48h10Tu:Leaf25  -0.12585    0.09569 252.00000  -1.315  0.18962   
## Treatment219Te48h1Tu:Leaf25   -0.22977    0.09569 252.00000  -2.401  0.01706 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 28 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
summary(mla_Tu)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(diffTu + 1) ~ Treatment2 * Leaf2 + (1 | Block2)
##    Data: new_total_box
## 
## REML criterion at convergence: -59.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3319 -0.3926 -0.0802  0.3748  3.1119 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Block2   (Intercept) 0.00000  0.0000  
##  Residual             0.03577  0.1891  
## Number of obs: 280, groups:  Block2, 2
## 
## Fixed effects:
##                                Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                   3.812e-02  5.981e-02  2.520e+02   0.637 0.524468
## Treatment210Te10Tu           -1.184e-01  8.458e-02  2.520e+02  -1.400 0.162797
## Treatment219Te1Tu            -7.500e-02  8.458e-02  2.520e+02  -0.887 0.376065
## Treatment21Te19Tu48h         -1.472e-01  8.458e-02  2.520e+02  -1.741 0.082929
## Treatment210Te10Tu48h        -1.154e-01  8.458e-02  2.520e+02  -1.364 0.173674
## Treatment210Te48h10Tu         4.201e-02  8.458e-02  2.520e+02   0.497 0.619878
## Treatment219Te48h1Tu         -3.390e-01  8.458e-02  2.520e+02  -4.008 8.06e-05
## Leaf23                        4.902e-04  8.458e-02  2.520e+02   0.006 0.995380
## Leaf24                       -2.032e-01  8.458e-02  2.520e+02  -2.402 0.017026
## Leaf25                       -4.833e-03  8.458e-02  2.520e+02  -0.057 0.954481
## Treatment210Te10Tu:Leaf23     1.201e-01  1.196e-01  2.520e+02   1.004 0.316407
## Treatment219Te1Tu:Leaf23     -5.645e-02  1.196e-01  2.520e+02  -0.472 0.637413
## Treatment21Te19Tu48h:Leaf23   1.064e-01  1.196e-01  2.520e+02   0.890 0.374532
## Treatment210Te10Tu48h:Leaf23  6.027e-02  1.196e-01  2.520e+02   0.504 0.614793
## Treatment210Te48h10Tu:Leaf23 -3.288e-02  1.196e-01  2.520e+02  -0.275 0.783627
## Treatment219Te48h1Tu:Leaf23   4.796e-01  1.196e-01  2.520e+02   4.010 8.02e-05
## Treatment210Te10Tu:Leaf24     1.025e-01  1.196e-01  2.520e+02   0.857 0.392153
## Treatment219Te1Tu:Leaf24      2.287e-01  1.196e-01  2.520e+02   1.912 0.057062
## Treatment21Te19Tu48h:Leaf24   3.930e-01  1.196e-01  2.520e+02   3.286 0.001162
## Treatment210Te10Tu48h:Leaf24  4.161e-01  1.196e-01  2.520e+02   3.479 0.000594
## Treatment210Te48h10Tu:Leaf24  1.753e-02  1.196e-01  2.520e+02   0.147 0.883606
## Treatment219Te48h1Tu:Leaf24   5.114e-01  1.196e-01  2.520e+02   4.276 2.71e-05
## Treatment210Te10Tu:Leaf25     2.398e-01  1.196e-01  2.520e+02   2.005 0.046033
## Treatment219Te1Tu:Leaf25      1.919e-02  1.196e-01  2.520e+02   0.160 0.872649
## Treatment21Te19Tu48h:Leaf25   1.129e-01  1.196e-01  2.520e+02   0.944 0.346168
## Treatment210Te10Tu48h:Leaf25 -7.036e-03  1.196e-01  2.520e+02  -0.059 0.953143
## Treatment210Te48h10Tu:Leaf25 -1.267e-01  1.196e-01  2.520e+02  -1.059 0.290719
## Treatment219Te48h1Tu:Leaf25   2.014e-01  1.196e-01  2.520e+02   1.684 0.093459
##                                 
## (Intercept)                     
## Treatment210Te10Tu              
## Treatment219Te1Tu               
## Treatment21Te19Tu48h         .  
## Treatment210Te10Tu48h           
## Treatment210Te48h10Tu           
## Treatment219Te48h1Tu         ***
## Leaf23                          
## Leaf24                       *  
## Leaf25                          
## Treatment210Te10Tu:Leaf23       
## Treatment219Te1Tu:Leaf23        
## Treatment21Te19Tu48h:Leaf23     
## Treatment210Te10Tu48h:Leaf23    
## Treatment210Te48h10Tu:Leaf23    
## Treatment219Te48h1Tu:Leaf23  ***
## Treatment210Te10Tu:Leaf24       
## Treatment219Te1Tu:Leaf24     .  
## Treatment21Te19Tu48h:Leaf24  ** 
## Treatment210Te10Tu48h:Leaf24 ***
## Treatment210Te48h10Tu:Leaf24    
## Treatment219Te48h1Tu:Leaf24  ***
## Treatment210Te10Tu:Leaf25    *  
## Treatment219Te1Tu:Leaf25        
## Treatment21Te19Tu48h:Leaf25     
## Treatment210Te10Tu48h:Leaf25    
## Treatment210Te48h10Tu:Leaf25    
## Treatment219Te48h1Tu:Leaf25  .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 28 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
Anova(mla_Te)
Anova(mla_Tu)
#### Contrasts
levels(as.factor(new_total_box$Treatment))
## [1] "10Te10Tu"    "10Te10Tu48h" "10Te48h10Tu" "19Te1Tu"     "19Te48h1Tu" 
## [6] "1Te19Tu"     "1Te19Tu48h"
mat_Leaf2<-matrix(ncol=10, nrow=7, 0)

#Compare Tu arriving first 10:10
mat_Leaf2[1,1]<-1
mat_Leaf2[2,1]<- -1


#Compare Te arriving first 10:10
mat_Leaf2[1,2]<-1
mat_Leaf2[3,2]<- -1

#Compare Te arriving first 19:1
mat_Leaf2[4,3]<-1
mat_Leaf2[5,3]<- -1

#Compare Tu arriving first 1:19
mat_Leaf2[6,4]<-1
mat_Leaf2[7,4]<- -1

#Compare Same arrival density effect Te
mat_Leaf2[1,5]<-1
mat_Leaf2[4,5]<- -1

#Compare Same arrival density effect Tu
mat_Leaf2[1,6]<-1
mat_Leaf2[6,6]<- -1

#Compare Tu arriving first (both densities)
mat_Leaf2[1,7]<-1
mat_Leaf2[2,7]<- -1
mat_Leaf2[6,7]<-1
mat_Leaf2[7,7]<- -1

#Compare Te arriving first (both densities)
mat_Leaf2[4,8]<-1
mat_Leaf2[5,8]<- -1
mat_Leaf2[1,8]<-1
mat_Leaf2[3,8]<- -1


#Compare 1:19Tu  (both timepoints)
mat_Leaf2[1,9]<-1
mat_Leaf2[2,9]<- 1
mat_Leaf2[6,9]<- -1
mat_Leaf2[7,9]<- -1

#Compare 19Te:1 (both timepoints)
mat_Leaf2[1,10]<-1
mat_Leaf2[3,10]<- 1
mat_Leaf2[4,10]<- -1
mat_Leaf2[5,10]<- -1

rownames(mat_Leaf2)<-c(levels(as.factor(df_wsingle$Treatment)))
colnames(mat_Leaf2)<-c("Tufirst_10:10", "Tefirst_10:10", "Tefirst_19:1", "Tufirst_1:19", "DensityTe19", "DensityTu19", "Tu arriving first", "Te arriving first", "Tu density", "Te density")
mat_Leaf2
##             Tufirst_10:10 Tefirst_10:10 Tefirst_19:1 Tufirst_1:19 DensityTe19
## 10Te10Tu                1             1            0            0           1
## 10Te10Tu48h            -1             0            0            0           0
## 10Te48h10Tu             0            -1            0            0           0
## 19Te1Tu                 0             0            1            0          -1
## 19Te48h1Tu              0             0           -1            0           0
## 1Te19Tu                 0             0            0            1           0
## 1Te19Tu48h              0             0            0           -1           0
##             DensityTu19 Tu arriving first Te arriving first Tu density
## 10Te10Tu              1                 1                 1          1
## 10Te10Tu48h           0                -1                 0          1
## 10Te48h10Tu           0                 0                -1          0
## 19Te1Tu               0                 0                 1          0
## 19Te48h1Tu            0                 0                -1          0
## 1Te19Tu              -1                 1                 0         -1
## 1Te19Tu48h            0                -1                 0         -1
##             Te density
## 10Te10Tu             1
## 10Te10Tu48h          0
## 10Te48h10Tu          1
## 19Te1Tu             -1
## 19Te48h1Tu          -1
## 1Te19Tu              0
## 1Te19Tu48h           0
test_la_Te<-testInteractions(mla_Te, fixed=c("Leaf2"), custom = list(Treatment2=mat_Leaf2), adjustment = "fdr" )

test_la_Tu<-testInteractions(mla_Tu, fixed=c("Leaf2"), custom = list(Treatment2=mat_Leaf2), adjustment = "fdr" )

test_la_Te
test_la_Tu

Calculating agreggation scores to test if changes between leaf, leaf age and direction

We will use C-score (package bipartite, which normalises c_score between 0 and 1) to test for aggregation between the two species for the different treatments

C score

For the C score we need a presence absence matrix per leaflet inside of a leaf

str(df)
## 'data.frame':    3214 obs. of  17 variables:
##  $ Box         : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Treatment   : chr  "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
##  $ Block       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Leaf        : int  2 2 2 2 2 2 3 3 3 3 ...
##  $ Leaflet     : int  1 2 3 1 2 3 1 2 3 4 ...
##  $ Direction   : Factor w/ 2 levels "Down","Up": 1 1 1 2 2 2 1 1 1 1 ...
##  $ Te          : int  4 2 1 6 6 3 0 0 4 0 ...
##  $ Tu          : int  1 2 5 1 12 6 1 2 0 2 ...
##  $ leaf_age    : chr  "old" "old" "old" "old" ...
##  $ Density     : chr  "10:10" "10:10" "10:10" "10:10" ...
##  $ Timing      : chr  "same time" "same time" "same time" "same time" ...
##  $ ratio       : num  4 1 0.2 6 0.5 0.5 1 0.5 4 0.5 ...
##  $ Total       : int  5 4 6 7 18 9 1 2 4 2 ...
##  $ proportionTe: num  0.8 0.5 0.167 0.857 0.333 ...
##  $ proportionTu: num  0.2 0.5 0.833 0.143 0.667 ...
##  $ Block2      : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Box2        : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
df$Leaf_leaflet<-sapply(c(1:length(df[,1])), function(x){paste(df$Leaf[x],df$Leaflet[x], sep="_")})
df$Leaf_leaflet_dir<-sapply(c(1:length(df[,1])), function(x){paste(df$Leaf[x],df$Leaflet[x], df$Direction[x], sep="_")})

### Summing everything per leaflet

aux_df2<-df %>%
  group_by(Treatment, Density, Timing , Block, Box, Leaf, Leaflet) %>%
     summarise(Te=sum(Te, na.rm=TRUE), Tu=sum(Tu, na.rm=TRUE)) %>%
  as.data.frame()
## `summarise()` has grouped output by 'Treatment', 'Density', 'Timing', 'Block', 'Box', 'Leaf'. You can override using the `.groups` argument.
aux_df2$Leaf_leaflet<-sapply(c(1:length(aux_df2[,1])), function(x){paste(aux_df2$Leaf[x],aux_df2$Leaflet[x], sep="_")})

aux_df2$Tu[which(aux_df2$Tu >1)]<-1
aux_df2$Te[which(aux_df2$Te >1)]<-1
str(aux_df2)
## 'data.frame':    1607 obs. of  10 variables:
##  $ Treatment   : chr  "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
##  $ Density     : chr  "10:10" "10:10" "10:10" "10:10" ...
##  $ Timing      : chr  "same time" "same time" "same time" "same time" ...
##  $ Block       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Box         : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Leaf        : int  2 2 2 3 3 3 3 3 4 4 ...
##  $ Leaflet     : int  1 2 3 1 2 3 4 5 1 2 ...
##  $ Te          : num  1 1 1 0 0 1 1 1 1 1 ...
##  $ Tu          : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ Leaf_leaflet: chr  "2_1" "2_2" "2_3" "3_1" ...
tu_df2<-pivot_wider(aux_df2[,c(1, 5,  9,10)], names_from=Leaf_leaflet, values_from=Tu)
te_df2<-pivot_wider(aux_df2[,c(1, 5, 8,10)], names_from=Leaf_leaflet, values_from=Te)

tu_df2$species<-"Tu"
te_df2$species<-"Te"

mat_02<-as.data.frame(rbind(tu_df2, te_df2))

mat_02
mat_02[is.na(mat_02)]<-0

c_score_mat<-expand.grid(Box=c(1:10), Treatment=levels(as.factor(df$Treatment)))
c_score_mat$cscore<-NA

str(mat_02)
## 'data.frame':    180 obs. of  21 variables:
##  $ Treatment: chr  "10Te10Tu" "10Te10Tu" "10Te10Tu" "10Te10Tu" ...
##  $ Box      : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ 2_1      : num  1 1 1 1 1 1 1 1 0 1 ...
##  $ 2_2      : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ 2_3      : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ 3_1      : num  1 0 1 0 0 0 1 1 1 1 ...
##  $ 3_2      : num  1 1 0 1 1 0 1 1 1 1 ...
##  $ 3_3      : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ 3_4      : num  1 1 1 1 1 1 1 0 1 1 ...
##  $ 3_5      : num  1 1 0 1 1 1 1 0 1 1 ...
##  $ 4_1      : num  1 1 1 1 1 1 1 1 1 0 ...
##  $ 4_2      : num  1 1 1 1 0 1 1 1 1 1 ...
##  $ 4_3      : num  1 1 1 1 1 1 1 1 0 1 ...
##  $ 4_4      : num  1 1 0 1 0 1 1 1 0 0 ...
##  $ 4_5      : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ 5_1      : num  1 0 1 1 1 1 1 1 1 1 ...
##  $ 5_2      : num  1 0 1 1 0 1 1 1 1 1 ...
##  $ 5_3      : num  1 0 1 1 1 1 1 1 1 1 ...
##  $ 5_4      : num  1 0 1 1 0 1 1 1 1 1 ...
##  $ 5_5      : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ species  : chr  "Tu" "Tu" "Tu" "Tu" ...
c_score_mat$cscore<-sapply(c(1:length(c_score_mat$Box)), function(x){
  C.score(t(subset(mat_02, Box==c_score_mat$Box[x] & as.character(Treatment)==as.character(c_score_mat$Treatment[x]))[,c(3:20)]), normalise = TRUE)
})


c_score_mat$Treatment2<-factor(c_score_mat$Treatment, levels=c("1Te19Tu", "10Te10Tu","19Te1Tu", "1Te19Tu48h","10Te10Tu48h","10Te48h10Tu","19Te48h1Tu"))

c_score_mat$Treatment3<-mapvalues(c_score_mat$Treatment2,c("1Te19Tu", "10Te10Tu","19Te1Tu", "1Te19Tu48h","10Te10Tu48h","10Te48h10Tu","19Te48h1Tu"), c("1:19\nsame time", "10:10\nsame time", "19:1\nsame time", "1:19\nTu first", "10:10\nTu first","10:10\nTe first", "19:1\nTe first") )

c_score_mat<-c_score_mat[-which(is.na(c_score_mat)),]

c_score_mat<-c_score_mat[-which(c_score_mat$cscore=="NaN"),]


c_score_mat<-droplevels(c_score_mat)

Figure S5

ggplot(c_score_mat, aes(x=Treatment3, y=cscore, fill=Treatment3))+
  geom_boxplot(na.rm=TRUE)+
  theme_ines+
  scale_fill_manual(values = c("#e0f3f8", "#91bfdb", "#4575b4",  "#fee090", "#ffffbf","#d73027", "#fc8d59"), name="Treatment")+
  scale_x_discrete(labels=c("","","","","","",""))+
  ylab("C-Score")+
  xlab("Treatment")+
  theme(axis.text.x = element_text(angle=60, vjust = 0.5))

save_plot("./Plots/C_score.pdf")

Test if C scores differ between treatments

shapiro.test(c_score_mat$cscore)
## 
##  Shapiro-Wilk normality test
## 
## data:  c_score_mat$cscore
## W = 0.62001, p-value = 4.254e-11
hist(c_score_mat$cscore)

descdist(c_score_mat$cscore+1, discrete = FALSE, boot=500)

## summary statistics
## ------
## min:  1   max:  2 
## median:  1 
## mean:  1.166668 
## estimated sd:  0.29382 
## estimated skewness:  2.014817 
## estimated kurtosis:  6.210139
plot(fitdist(c_score_mat$cscore+1, "gamma"))
## $start.arg
## $start.arg$shape
## [1] 16.03821
## 
## $start.arg$rate
## [1] 13.74702
## 
## 
## $fix.arg
## NULL

plot(fitdist(c_score_mat$cscore, "beta", "mme"))
## $start.arg
## $start.arg$shape1
## [1] 0.1060946
## 
## $start.arg$shape2
## [1] 0.5304664
## 
## 
## $fix.arg
## NULL

plot(fitdist(c_score_mat$cscore+1, "lnorm"))
## $start.arg
## $start.arg$meanlog
## [1] 0.1294786
## 
## $start.arg$sdlog
## [1] 0.2097881
## 
## 
## $fix.arg
## NULL

fit_dist_cscore<-summary(fitdist(c_score_mat$cscore, "beta", "mme"))
## $start.arg
## $start.arg$shape1
## [1] 0.1060946
## 
## $start.arg$shape2
## [1] 0.5304664
## 
## 
## $fix.arg
## NULL
str(fit_dist_cscore)
## List of 20
##  $ estimate   : Named num [1:2] 0.106 0.53
##   ..- attr(*, "names")= chr [1:2] "shape1" "shape2"
##  $ method     : chr "mme"
##  $ sd         : NULL
##  $ cor        : NULL
##  $ vcov       : NULL
##  $ loglik     : num Inf
##  $ aic        : num -Inf
##  $ bic        : num -Inf
##  $ n          : int 59
##  $ data       : num [1:59] 1 0 0 0.333 0.444 ...
##  $ distname   : chr "beta"
##  $ fix.arg    : NULL
##  $ fix.arg.fun: NULL
##  $ dots       : NULL
##  $ convergence: num 0
##  $ discrete   : logi FALSE
##  $ weights    : NULL
##  $ ddistname  : chr "dbeta"
##  $ pdistname  : chr "pbeta"
##  $ qdistname  : chr "qbeta"
##  - attr(*, "class")= chr [1:2] "summary.fitdist" "fitdist"
c_score_mat$Block<-mapvalues(c_score_mat$Box, c(1,2,3,4,5,6,7,8,9,10), c(1,1,1,1,1,2,2,2,2,2))
c_score_mat$Block<-as.factor(c_score_mat$Block)

cs2<-glm((cscore+1)~Treatment, family=Gamma, data=c_score_mat)

summary(cs2)
## 
## Call:
## glm(formula = (cscore + 1) ~ Treatment, family = Gamma, data = c_score_mat)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.35779  -0.13385  -0.03956   0.00000   0.56236  
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           0.68354    0.06092  11.220 1.68e-15 ***
## Treatment10Te10Tu48h  0.14278    0.08560   1.668  0.10135    
## Treatment10Te48h10Tu  0.10912    0.08390   1.301  0.19915    
## Treatment19Te1Tu      0.31646    0.09834   3.218  0.00222 ** 
## Treatment19Te48h1Tu   0.25675    0.09865   2.603  0.01203 *  
## Treatment1Te19Tu      0.15264    0.08393   1.819  0.07473 .  
## Treatment1Te19Tu48h   0.27742    0.09007   3.080  0.00331 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.04766633)
## 
##     Null deviance: 2.9115  on 58  degrees of freedom
## Residual deviance: 2.1399  on 52  degrees of freedom
## AIC: 3.3808
## 
## Number of Fisher Scoring iterations: 4
Anova(cs2, type=2)
testInteractions(cs2, pairwise="Treatment", adjustment = "fdr")

Order of arrival affects probability of coexistence?

Calculating FD from Godoy and Levine 2014

The demographic parameter in our case is lambdai -1 / lambdaj -1 and the competitive response the alphas which are arranged in a different way to calculate niche diff.

aij describes the per capita effect of species j on species i –> so in the script, in the results the row is the focal so the Te inter is the effect that Tu has on Te, and Tu inter vice versa.

fitness ratio = ((nj-1)/(ni-1))* sqrt(aij/ajj * aii/aji) Niche overlap = sqrt(aij/ajj* aji/aii)

Don’t forget that niche differences is 1- niche overlap

The data is generated by the script alpha_estimate_fixed_final.Rmd, so please run this script before hand

df_alpha_lambda<-read.csv("./Estimate_Lambda_fixed_allTreats/Results_all.csv", header=TRUE)

ggplot(df_alpha_lambda, aes(x=Comp, y=Alpha, colour=Arrival, fill=Arrival))+
  facet_grid(.~Species)+
  geom_errorbar(aes(ymin=Alpha_Lower, ymax=Alpha_Upper), position=position_dodge(0.5), size=0.75)+
  geom_point(size=2, position=position_dodge(0.5), shape=21, colour="black")+
  geom_line(aes(group=Arrival), position=position_dodge2(0.5))+
  scale_colour_manual(values = c("#4575b4","#d73027", "#fee090"))+
  scale_fill_manual(values =c("#4575b4","#d73027", "#fee090"))+
  theme_ines+
  ylab("Competitive ability")+
  xlab("")+
  theme(legend.position = "bottom")

save_plot("./Estimate_Lambda_fixed_allTreats/alpha.pdf")

ggplot(df_alpha_lambda, aes(x=Species, y=Lambda, colour=Arrival, fill=Arrival))+
   geom_hline(yintercept =2, colour="lightgray")+
  geom_hline(yintercept =3.5, colour="lightgray")+
  geom_errorbar(aes(ymin=Lambda_Lower, ymax=Lambda_Upper), position=position_dodge(0.5), size=0.75)+
  geom_point(size=2, position=position_dodge(0.5), shape=21, colour="black")+
  scale_colour_manual(values = c("#4575b4","#d73027", "#fee090"))+
  scale_fill_manual(values = c("#4575b4","#d73027", "#fee090"))+
  theme_ines+
  ylab("Lambda Growth Rate")+
  xlab("")+
  theme(legend.position = "bottom")

save_plot("./Estimate_Lambda_fixed_allTreats/lambda.pdf")


coex_df2<-as.data.frame(matrix(nrow=3, ncol=7))
colnames(coex_df2)<-c("Treatment", "NO", "NO_L", "NO_U", "FD", "FD_L","FD_U")

coex_df2$Treatment<-c("same time","Te first","Tu first")
str(coex_df2)
## 'data.frame':    3 obs. of  7 variables:
##  $ Treatment: chr  "same time" "Te first" "Tu first"
##  $ NO       : logi  NA NA NA
##  $ NO_L     : logi  NA NA NA
##  $ NO_U     : logi  NA NA NA
##  $ FD       : logi  NA NA NA
##  $ FD_L     : logi  NA NA NA
##  $ FD_U     : logi  NA NA NA
str(df_alpha_lambda)
## 'data.frame':    12 obs. of  9 variables:
##  $ Species     : chr  "Te" "Tu" "Te" "Tu" ...
##  $ Arrival     : chr  "same time" "same time" "Te first" "Te first" ...
##  $ Comp        : chr  "Inter" "Inter" "Inter" "Inter" ...
##  $ Lambda      : num  3.6 1.73 3.6 1.73 3.6 ...
##  $ Lambda_Upper: num  4.18 1.9 4.18 1.9 4.18 ...
##  $ Lambda_Lower: num  3.02 1.56 3.02 1.56 3.02 ...
##  $ Alpha       : num  0.0532 0.0261 0.0586 0.0455 0.0679 ...
##  $ Alpha_Upper : num  0.0648 0.0368 0.0693 0.0731 0.0767 ...
##  $ Alpha_Lower : num  0.0417 0.0153 0.048 0.0178 0.0591 ...
head(df_alpha_lambda)
### Important!
#The order here is Te inter (aji), Tu inter (aij), Te intra (ajj), Tu intra (aii) for each treatment

coex_df2$NO<-sapply(c(1:3), function(x){
  aux<-subset(df_alpha_lambda, as.character(Arrival)==as.character(coex_df2$Treatment[x]))
  
  #Assuming that the order is always te tu, te tu, inter inter, intra intra
  nd<-(sqrt(aux$Alpha[2]/aux$Alpha[3]* aux$Alpha[1]/aux$Alpha[4]))
  
  nd
})

coex_df2$NO_L<-sapply(c(1:3), function(x){
  aux<-subset(df_alpha_lambda, as.character(Arrival)==as.character(coex_df2$Treatment[x]))
  
  #Assuming that the order is always te tu, te tu, inter inter, intra intra
  nd<- sqrt(aux$Alpha_Lower[2]/aux$Alpha_Lower[3]* aux$Alpha_Lower[1]/aux$Alpha_Lower[4])
  
  nd
})

coex_df2$NO_U<-sapply(c(1:3), function(x){
  aux<-subset(df_alpha_lambda, as.character(Arrival)==as.character(coex_df2$Treatment[x]))
  
  #Assuming that the order is always te tu, te tu, inter inter, intra intra
  nd<-(sqrt(aux$Alpha_Upper[2]/aux$Alpha_Upper[3]* aux$Alpha_Upper[1]/aux$Alpha_Upper[4]))
  
  nd
})

#fd= ((nj-1)/(ni-1))* sqrt(aij/ajj * aii/aji)

coex_df2$FD<-sapply(c(1:3), function(x){
  aux<-subset(df_alpha_lambda, as.character(Arrival)==as.character(coex_df2$Treatment[x]))
  
  #The order here is Te inter (aji), Tu inter (aij), Te intra (ajj), Tu intra (aii) 
  fd<-((aux$Lambda[1]-1)/(aux$Lambda[2]-1))*(sqrt(aux$Alpha[2]/aux$Alpha[3])* sqrt(aux$Alpha[4]/aux$Alpha[1]))
  
  fd
})

coex_df2$FD_L<-sapply(c(1:3), function(x){
  aux<-subset(df_alpha_lambda, as.character(Arrival)==as.character(coex_df2$Treatment[x]))
  
  #Assuming that the order is always te tu, te tu, inter inter, intra intra
  fd<-((aux$Lambda_Lower[1]-1)/(aux$Lambda_Lower[2]-1))*(sqrt(aux$Alpha_Lower[2]/aux$Alpha_Lower[3])* sqrt(aux$Alpha_Lower[4]/aux$Alpha_Lower[1]))
  
  fd
})

coex_df2$FD_U<-sapply(c(1:3), function(x){
  aux<-subset(df_alpha_lambda, as.character(Arrival)==as.character(coex_df2$Treatment[x]))
  
  #Assuming that the order is always te tu, te tu, inter inter, intra intra
  fd<-((aux$Lambda_Upper[1]-1)/(aux$Lambda_Upper[2]-1))*(sqrt(aux$Alpha_Upper[2]/aux$Alpha_Upper[3])* sqrt(aux$Alpha_Upper[4]/aux$Alpha_Upper[1]))
  
  fd
})


coex_df2$ND<-1-coex_df2$NO
coex_df2$ND_L<-1-coex_df2$NO_L
coex_df2$ND_U<-1-coex_df2$NO_U

bound_df<-data.frame(niche_overlap=c(seq(0,3, 0.05))) # creating a vector with niche overlap
bound_df$niche_diff<-(1-bound_df$niche_overlap) # calculating stabilizating differences from niche overlap 1-rho
bound_df$fitness_differences_sp_1<-(1/bound_df$niche_overlap) # solid line in your graph this is ok
bound_df$fitness_differences_sp_temp<- 1-bound_df$fitness_differences_sp_1 #this is an intermediate step to see the differences above one 
#which is later incorporated into the 2 species
bound_df$fitness_differences_sp_2<- 1+ bound_df$fitness_differences_sp_temp# dashed line in your graph which is NOT OK in your original graph
#Here I added the difference calculated in the previous step to make the curves simetric. they must be simmetric!
#otherwise one species has more chance of priority effects than the other.
bound_df<-bound_df[, -4]
#remove the intermediate step 
str(bound_df)
## 'data.frame':    61 obs. of  4 variables:
##  $ niche_overlap           : num  0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 ...
##  $ niche_diff              : num  1 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 ...
##  $ fitness_differences_sp_1: num  Inf 20 10 6.67 5 ...
##  $ fitness_differences_sp_2: num  -Inf -18 -8 -4.67 -3 ...
colnames(bound_df)<-c("NO", "ND", "FD", "FD2")

Fig 2 - Coexistence plot

ggplot(coex_df2,aes(x=ND, y=FD, colour=Treatment) )+
  geom_hline(yintercept = 1 ,colour="lightgray")+
  geom_vline(xintercept = 0 , colour="lightgray")+
  geom_ribbon(data=bound_df, aes(x=ND, ymin=FD,ymax=FD2), colour="Black", fill="lightgrey", alpha=0.5)+
  geom_line(data=bound_df, aes(x=ND, y=FD), colour="black")+
  geom_line(data=bound_df, aes(x=ND, y=FD2),colour="black")+
  geom_errorbar(data=coex_df2, aes(ymin=FD_L, ymax=FD_U), colour="black", width=0.1)+
  geom_errorbarh(data=coex_df2, aes(xmin=ND_L, xmax=ND_U), colour="black", height=0.1)+
  geom_point(data=coex_df2, size=2.75, shape=21, aes(fill=Treatment), colour="black")+
  ylab(c(""))+
  xlab(c(""))+
  #scale_x_discrete(labels=c())+
  scale_colour_manual(values = c("#4575b4","#d73027", "#fee090"), name="")+
  scale_fill_manual(values = c("#4575b4","#d73027", "#fee090"), name="")+
  scale_x_continuous(breaks = seq(-1,1, 0.2), expand=expansion(c(0,0)), limits = c(-1.5, 1))+
  ylim(c(-1.5,5))+
  theme_bw()+
  theme_ines+
  theme(legend.position = c(0.9,0.9), legend.background = element_rect(fill=NA), legend.text = element_text(size=10), axis.text = element_text(colour="black"))+
  geom_richtext(label="<span style='color:black'>Priority effects", x=-0.82, y=1.1, size=4, color=NA, fill=NA)+
  geom_richtext(label="<span style='color:black'>*T. urticae* excluded", x=0, y=2.5, size=4, color=NA, fill=NA)+
  geom_richtext(label="<span style='color:black'>*T. evansi* excluded", x=0, y=0, size=4, color=NA, fill=NA)+
  geom_richtext(label="<span style='color:black'>Coexistence", x=0.56, y=1.1, size=4, color=NA, fill=NA)+
  expand_limits(x = c(-1, 0.75), y = 0)+
  coord_cartesian(ylim = c(0, 5), xlim = c(-1,0.7))
## Warning: Removed 10 row(s) containing missing values (geom_path).

## Warning: Removed 10 row(s) containing missing values (geom_path).

save_plot("./Estimate_Lambda_fixed_allTreats/coex.pdf")
## Warning: Removed 10 row(s) containing missing values (geom_path).

## Warning: Removed 10 row(s) containing missing values (geom_path).